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ABSTRACT 

Compelling evidence associates the nuclei of active galaxies and massive star- 
bursts. The symbiosis between a compact nuclear starburst stellar cluster and a mas- 
sive black hole can self-consistcntly explain the properties of active nuclei. The young 
stellar cluster has a profound effect on the most important observable properties of 
active galaxies through its gravity, and by mass injection through stellar winds, super- 
novae and stellar collisions. This mass-loss, injected throughout the nucleus, creates 
a hot nuclear interstellar medium (nISM). The cluster both acts as an optically-thin 
fuel reservoir and enriches the nISM with the products of nucleosynthesis. The nISM 
flows under gravitational and radiative forces until it leaves the nucleus or is accreted 
onto the black hole or accretion disc. 

The radiative force exerted by the black hole-accretion disc radiation field is 
not spherically symmetric. This results in complex flows in which regions of inflow 
can coexist with high Mach number outflowing winds and hydrodynamic jets. We 
present two-dimensional hydrodynamic models of such nISM flows, which are highly 
complex and time variable. Shocked shells, jets and explosive bubbles are produced, 
with bipolar winds driving out from the nucleus. Our results graphically illustrate why 
broad emission line studies have consistently failed to identify any simple, global flow 
geometry. The real structure of the flows is inevitably yet more complex. 

The structure of these nISM flows is principally determined by two dimensionless 
quantities. The first is the magnitude of the stellar cluster velocity dispersion relative 
to the sound speed in the nISM. These speeds measure the gravitational and thermal 
energies in the nlSM respectively, and, therefore, whether the gas is initially bound, 
or escapes in a thermal wind. The second parameter is the Mach number of the ill- 
collimated nISM flow which is driven away from the central black hole. We discuss a 
two-parameter classification based on this observation which, in future papers, we will 
relate to empirical classifications. 

The interplay between the nucleus and the wider galaxy depends critically on 
the exchange of radiative and mechanical energy. The outbound mechanical energy 
transfer is governed by the nuclear stellar cluster. Active galactic nuclei will only 
be understood once the symbiotic relationships between the black hole, the stellar 
cluster, and the galaxy are considered, ft is impossible to treat correctly any isolated 
component. Our conceptually simple and self-consistent Symbiotic Model explains the 
observed complexity of active galaxies without ad hoc measures. 
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There is strong evidence for black holes with a wide range in 
mass in the centres of galaxies, and unambiguous evidence 
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for dense stellar clusters in galactic cores. Activity in the nu- 
clei of galaxies appears to be intimately linked not only to 
these supermassive black holes (BH) and nuclear starbursts, 
but also to interactions between galaxies. The evidence for 
this link - both direct and circumstantial - has been re- 



viewed in detail by Perry (1993a, 1999; see also Pilippenko 
1992()7Jhe exact nature of the link remains uncertain: inter- 
actions, starbursts and activity certainly can not be equated, 
and not all systems with these components are active. In 
partic ular, galactic interactions do not always result in ac- 
tivity ( Sanders fc Mirabel 1996 ) . However, the evidence does 
unequivocally demonstrate that the active nucleus is not a 
'parasite' upon a 'host' galaxy, but rather that activity is a 
galaxy- wide phenomenon which is most clearly identified by 
the dramatic events in the galaxy nucleus. In this paper, we 
consider the hydrodynamics of this basic nuclear system. 

The most dramatic of the nuclear phenomena is the 
highly luminous (lO'^" to lO'^^L©), time- variable, broad- 
band (radio to 7-ray) radiation emitted from regions as 
small as a few light-seconds. This radiation is accompanied 
in many AGN by broad emission lines from heavy elements 
in cold (^ IQ-'K) dense {'^ 10^ to 10^^ cm'^*) gas. In the 
unified model of active nuclei, the broad line region is a 
constituent of all active nuclei, although it is not observed 
in some as a result of continuum beaming or inter vening 
obscuration ( Krolik fc Begelman 1986 ; Barthel 1989 ). Line 
widths, which reflect local gas velocities, are typically sev- 
eral lO'^kms^^. Average equivalent widths are remarkably 
consistent across orders of magnitude differences in total 
bolometric luminosity (although with a marked scatter be- 
tween individual objects at a given luminosity) , indicating a 
tight relationship between the production of cold, condensed 
radiating gas structures and the production of the central 
luminosity itself. 

Although accretion onto black holes seems clearly to 
dominate the physics of the innermost regions of AGN 
(within ~ 100 Schwarzschild radii), and to provide a sat- 
isfying explanation of the radiation processes which create 
the luminous broad-band continuum, pure black hole accre- 
tion models have always suffered from an inability to explain 
the structure of the line emitting regions and had great dif- 
ficulty in explaining the fuelling, before the importance of 
stars in the wider parsec sized core was recognised. 

The stars which play the dominant role in the energet- 
ics and mass budget of galactic nuclei are young, mass-losing 
stars: these are most abundant and important in just those 
nuclear stellar clusters - starbursts - which are thought to 
be involved in activity. Current theoretical models linking 
activity and starbursts come in two vari ants. In this paper, 
we explore the structure of the model ( |^erry 1992 ) b ased 
on the early work of Perry & Dyson (1985, hereafter PD) 
which unites starbursts and black holes as compact symbi- 
otic systems in the nuclei of active galaxies (CoUin-Souffrin 
et al 1988, hereafter CDMP; [Norman fc Scoville 1988| ; [Perry 
1992; V/illiams & Perry 1994). In contrast, the pure star- 



burst model equates the two phenomena, consider ing AGN 
to b e but extreme example s of nuclear starbursts ( Perlevich 



199C 



C^id Fernandes 1997). Interestingly, radiative shocks 



ergetics of starbursts, understanding the dynamic link be- 
tween black holes and starbursts requires careful analysis of 
the stellar- and hydro-dynamics of the combined system. In 
this paper we present the first comprehensive treatment of 
the hydrodynamics of such a symbiotic system. 

In the rest of this introduction, we outline the physical 
context and assumptions of our model. These assumptions 
will be discussed in more detail in subsequent sections. 

1.1 Symbiotic models 

In the Symbiotic Model of AGN, the nuclear starburst stellar 
cluster feeds the black hole through the nuclear interstellar 
medium (nISM) created by stellar mass-loss. The formation 
of shocks within the nISM is an inescapable consequence 
of the interaction between individual (mass-loading) stars 
and the global nISM. It is important to emphasise that the 
presence of shocks in the nISM is a model independent phe- 
nomenon and is a direct consequence of basic thermodynam- 
ics, stellar dynamics and hydrodynamics in galactic nuclei. 
In QSOs, the shocks around individual supernovae are suf- 
ficiently large that the shocked nISM gas has time to cool 
radiatively into thermal equilibrium with the central radia- 
tion field at ~ 10* K, so forming broad line emitting struc- 
tures (PD). The thermal pressure in such shock confined 
structures equals the stagnation pressure of the flow, which 
agrees well with the pressures deduced from observations of 
the emission lines. PD showed that using the observed equiv- 
alent widths of the high-ionization lines (HILs) to deduce the 
average properties of the nISM, the implied flow rates in the 
nISM are sufficient to fuel the black hole accretion. 

This satisfyingly self-consistent picture did not initially 
explain the two-(or more)-component nature of the broad 
emission line region (BELR): the low-ionization lines (LILs) 
are observed to be both red-shifted with respect to the 
HILs, and to require a harder ionizing continuum. However, 
CDMP then showed that the nISM fiow can back-scatter X- 
rays which provide 'deep heating' of the accretion disc suf- 
ficient to produce the LILs. This two-component model re- 
solved several outstanding questions concerning the BELR. 
The transient nature of the BELR clouds solved the 'cloud 
confinement problem'. The existence of two distinct regions 
of broad line emission also gave a simple explanation for the 
observed systematic velocity shifts and multiple component 
structure of the broad line profiles (Corbin fc Francis 1994 ; 
Brotherton et al. 1994; Arav et al. 1998). Our results, in this 



paper, add another dimension to this understanding of the 
structure of the BELR: we find that the nISM itself forms 
distinct multi-component regions of differing characteristic 
excitation levels and velocities. 

The structure of the emission line region which emerges 
is that of distinct local clumps of line-emitting gas confined 
by transient shocks. Recently, Arav et al. used the new gen- 
eration of observations to update the limits which smooth 
line profiles put on the number of such centres of emission 
(e.g. Atwood, Baldwin fc Carswell 1982). Their analysis was 



around supernovae play a central role in both theories. We 
shall refer to the starburst-black hole paradigm for AGN as 
the Symbiotic Model. 

Since powerful winds and supernovae dominate the en- 



based on the fact that discrete emitting units give a smooth 
line profile only when the number of emitting units is very 
large. Their method is sensitive to structures with veloc- 
ity widths between 50 and 200kms~^; within this range 
they find essentially no perturbations in the wings of the 
Balmer lines, to very tight limits. This rules out many of 
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the simplest scenarios for the formation of the lines by dis- 
crete clumps. Arav et al. were able to rule out emission 
from less than 3 x l(f smooth shell sources with expan- 
sion velocity lOOkms"^. However, such sources, with small 
velocity width and perfect symmetry, are not representa- 
tive of shock-confined emission line clouds in galactic nuclei. 
Structures which are intrinsically smoother in emission than 
those they specifically considered, such as accretion discs 
(which are thought on other grounds to be important con- 
tributors to the Balmer emission, CDMP), or strong wind 
sources and bowshocks {v ^ 200 kms~^) are not ruled out by 
these observations. The sonic velocity in the nISM is at least 
400kms~^; this is the minimum expected shock velocity in 
the region. In fact, characteristic velocities of the shocks 
which arise in symbiotic nuclei are xlO^kms"^ [the combi- 
nation of stellar ejection velocities (often 3 x 10'' kms~^) 
and the stellar orbital velocities (^ 600 kms^^ within the 
central parsec)]. Furthermore, the shocks are non-uniform: 
they are either open bow shocks with streaming tails, or 
closed shells distorted by their interaction with the wind. 
The density, radial velocity and surface emissivity - for any 
shock - all vary significantly around the shock- front, gener- 
ating broad, complex lines with wings which are much less 
steep than those Arav et al. assumed. As Arav et al. pointed 
out, any flow which produces smoothly varying structures on 
velocity scales of 500 - 1000 km s~^ is allowed by the obser- 
vations. We therefore conclude that their elegant analysis 
provides support for the PD model, and underlines the im- 
portance of hydrodynamic models of the BELR. 

1.2 Physical Ingredients 

The physical parameters characterising the centre of a 
galaxy are the mass of the black hole, and the mass, stellar 
density proflle, IMF and age of the stellar cluster. The IMF 
and age govern the mass-loss rate from stars due to stellar 
evolution; collisions and tidal disruption are negligible ex- 
cept in the central regions of the densest clusters. The mass 
lost from the stars inevitably creates a nuclear interstellar 
medium (nISM). If there is a black hole accretion flow, either 
from this nISM or from some other reservoir (e.g. flow into 
the nucleus from a galactic bar), then radiation processes 
near the black hole can generate a powerful broad-band con- 
tinuum (althou gh not if the flow int o the black hole is hot or 
optic ally thick (Narayan & Yi 1994; Blandford & Begelman 



199!:)). For the present, we take this as an additional input 
to our simulations. The system which we consider in this 
paper is one in which such accretion processes are effective. 

The stellar cluster surrounding the black hole helps to 
determine the structure of the nISM flow on the scale of the 
BELR, through gravity and distributed stellar mass loss. 
The strength of this mass-loading is principally determined 
by the age and initial mass function (IMF) of the stellar clus- 
ter. We shall assume that the cluster is young and dominated 
by high mass stars, and take our basic paramet ers fro m the 
stellar synthesis models of Williams and Perry (1994). 



Mass input from the stellar cluster is discrete and local. 
Because we are interested in the generic properties of the 
global flows, we will concentrate on the luminous QSOs. In 
such models, the mass processing rates are high enough that 
we can average over the numerous small scale structures re- 
sulting from the interaction between individual stars and 



supernovae with the global flow. Future papers will investi- 
gate the detailed structure of individual mass input sources, 
and the effects of the discrete mass input on the global flow. 

We have focussed on mass injection from stars within 
the cluster rather than from the external ISM, primarily for 
two reasons. Importantly, the mass flux from the external 
ISM is unlikely to be able to provide a significant fraction 
of the fuelling require ments of the AGN as a continuous 
source (Shlosman et al. 1990). The wind from the nucleus is 



in some cases strong enough to drive the external ISM away. 
Also, this is the situation described in PD, where the shocks 
around the mass-loading centres also generate the cool gas 
which emits the observed emission lines. Therefore, rather 
than include the additional free parameters specifying the 
distribution of the ambient ISM at the outer boundary of our 
models in this first paper, we have chosen to concentrate on 
the case where the phase of accretion from tens or hundreds 
of parsecs to within a parsec of the central black hole has 
already ceased. 

We have assumed that any relativistic jet driven from 
close to the central black hole is absent, or at least narrow, 
transferring little momentum to the surrounding global flow. 
If this assumption is reasonable, our models apply equally 
to radio-loud and radio-quiet objects. Radio structures on 
parsec scales are seen in radio-quiet objects as well as radio- 
loud ones: the explicit inclusion of jets in our models would 
be an interesting extension. 

1.3 Formation scenarios 

The starburst-black hole system is in a symbiotic relation- 
ship with the wider galaxy. We study the structure of the 
nuclear system in its active phase in this paper. We leave the 
vexed questions of the formation of the nuclear system and 
the details of its symbiosis with the wider galaxy to later 
work. This follows the tradition well established in stellar 
physics: first understand structure, then evolution, lastly for- 
mation. Nonetheless, there are clear indications, both from 
obse rvations an d from theory, of possible formation scenar- 
ios (Perry 1999). Interactions between galaxies disturb the 
equilibrium of the ISM of the galaxies, and gas can flow to- 
wards the nucleus where it can form a nuclear starburst. The 
size scale of the starburst depends on whether the gas can 
shed its angular mome ntum, and o n whether it hangs up at 



a Lindblad resonance (Lamb 1999) 



Direct optical imaging of galaxies with bright active nu- 



ing the morphological types (Punlop ct al. 1993; McLeod & 


Rieke 1994a|; 1994b 


; Taylor et al. 1996|; |Bahcall, Kirhakos 


& Schneider 1997; 


Boyce 1997; McLeod 1997). Radio and 



quasar galaxies appear to be large and luminous, and to 
show signs of past merger events. Galaxies with radio-loud 
active nuclei are essentially all spheroidal, whereas those 
with radio-quiet cores are approximately equally divided 
between spheroidal and disc-dominated structures (Kukula 



et al. 1997a). This may indicate e.g. different interaction his- 
tories for radio-loud and radio-quiet systems. 

There is much circumstantial evidence that similar en- 
vironmental factors influence the frequency of both active 
nuclei and starburst galaxies (Neff & Hutchings 1992). Cool- 
ing flows around massive elliptical galaxies may provide a 
source of fuel for both phenomena ( Binney fc Tabor 199f: ; 
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Ciotti fc Ostriker 19971; |Bremer, Fabian fc Crawford 1997t 



Friaga fc Terlevich 1998[). Ultraluminous IR galaxies may 



represent a transition phase as dense shrouds of dust lift to 
reveal bright optical quasar nuclei, since the environments, 
molecular gas content, near-IR colours and optical emission 



line diagnostics of the two classes are similar ( Banders et al 
198S)- Processed material with near-solar metallicities in the 



nISM is seen in AGN even at the highest redshifts. At red- 
shifts 2 < 2 < 3, absorption line studies reveal that most 
of the gas in g a laxies typically has Z/7iq « 0.01-0.1 (Pet- 
tini et ll. 1994 



1997 ). To reach the high metal abundances 



observed in QSOs after short cosmic times, massive star for- 
mation, probably m situ, m ust be intimatelv linked to the 



evolution of active galaxies ( Perry 1992 ; Ferland et al. 1996 
Korista et al. 1996). Populations of star-forming gala xies at 



2 ~ 3 have now been discovered (Steidel et al. 199C). This 
redshift range is hence a critical period in t he evolution o f 
both normal and quasar galaxy populations (McLeod 1997). 

Whatever the likely formation scenario, the nISM in our 
models reaches a structure characteristic of its long term 
behaviour - equilibrium, oscillations or explosions - on dy- 
namical timescales typically a few thousand years. These 



timescales arc m uch shorter than any plau sible formation or 



evolution times (Williams & Perry 1994), and so we treat 



the structure of the system independently from its formation 
and evolution. 



1.4 Hydrodynamics of the nISM 

Our aim in this paper is to explore the hydrodynamics of 
those nuclear bl ack-ho le /stellar cluster systems which are 
active. Williams (1999) summarises previous hydrodynamic 
models of the flow structure of active nuclei. Our study is 



unique in that it treats the distributed mass input through 
the nucleus and the gravitational field of the stellar cluster. 

We treat the global hydrodynamics of the nISM, as it 
interacts with the central black hole, the stellar cluster and 
the accretion disc. These physical components are essential 
to the starburst-black hole model of AGN - and to most 
other dynamical models of active nuclei. We have chosen 
to simplify the properties of this system as far as possible. 
This approach allows us focus on the most important proper- 
ties of the global hydrodynamic interactions between stellar 
cluster, black hole and accretion disc. 

Independent of the detailed model, it is possible to es- 
timate characteristic values of the flow density and velocity 
required for the accreting nISM flow to be able to fuel the 
AGN. PD used such simple scaling arguments to discuss the 
structure of the nISM flow, and to investigate the conse- 
quences for the fuelling of AGN and for the creation of the 
BELR. If the local nISM is too dense, it obscures the sys- 
tem; if it is too diffuse, it fails to provide the necessary fuel 
for the AGN. Those flows consistent with these constraints 
can be shown to have ram pressures similar to the pressures 
deduced from observations of the broad line clouds. The ram 
pressure of the nISM will therefore have an important dy- 
namical effect on the BELR clouds, whatever their confine- 
ment mechanism. Although these simple arguments provide 
interesting overall constraints, many important properties 
of the AGN depend on the structure of the flow, which PD 
did not address. Detailed hydrodynamic modelling is nec- 
essary to determine these flow structures. More broadly, it 



is clearly necessary to consider the hydrodynamics of the 
nISM in order to develop fully any model of AGN. 

In the absence of a black hole, in a dense stellar cluster 
with negligible external pressure, gas accretes onto a hydro- 
static core, and a weak wind forms (cf. Section]^). The core 
grows until it reaches the sonic radius, and the wind from 
the cluster increases to an equilibrium in which it trans- 
ports the mass injection from the stellar cluster out into the 
wider galaxy. Once a black hole is added to the system, with 
an accretion-generated radiation field, many more complex 
behaviours can occur, some of which we investigate in this 
paper. Broadly, however, the black hole can accrete matter 
at the centre of the fiow, weakening the eventual wind, and 
the radiation field may drive it away. In the absence of radia- 
tive driving, the sonic (or Parker) radius is GM/2cl — 1.4 pc 
for a 10^ Mq BH mass and gas at a temperature of 10*^ K, 
which is at least as large as the broad emission line region 
even in QSOs. The velocities in pure thermally driven winds 
only reach ~ 360kms~^ (the sound speed in the nISM) at 
such radii. Such low velocities are in confiict with the strong 
evidence for outfiowing gas in and around AGN at far higher 



velocities, most clearly in bro ad absorption line QSOs ( Turn- _ 
shek 1988, psterbrock 1991 ). If radiative driving occurs at 



all angles, then the hole may well be starved of fuel. This 
will reduce the luminosity, giving rise to a limiting feedback 
process. However, there is strong evidence of beaming in all 
AGN. This asymmetry in the radiation field means that the 
radiation force can be greater than the gravitational attrac- 
tion of the black hole on the axis of the accretion disc, but 
not in its plane. In such cases, mass can fall inwards towards 
the plane and be driven away on axis. Such flows are typical 
of our simulations. 

The structures and the global kinematics of these nISM 
flows are likely to be generic features. However, we will see 
that small-scale features such as gas cooling may play an im- 
portant role in determining the observable properties of the 
flow. Where cooling is significant, the average opacity of the 
gas will be increased, which will re-scale the importance of 
the central radiation field. The central broad-band (radio to 
X- and 7-rays) radiation field is powered by mass accretion 
through a disc and down into the black hole. We model this 
radiation field as the sum of an unresolved point source and 
an accretion disc component which follows a simple cosine 
law. We concentrate on models appropriate to QSOs, since 
there is wide agreement that they accrete at or near the 
Eddington limit. In contrast, the lower luminosity Seyferts 
appear to span a much wider range in mass-to-light ratio. 
We discuss the consequences of scaling our models to Seyfert 
properties in Section 6.5. 

Because we do not model the accretion disc explicitly, 
the effective Eddington ratio, /Edd, is an arbitrary input pa- 
rameter. A goal of future research is to determine the rela- 
tionship between /Edd and the properties of the nISM, since 
this must be self-consistent (i.e., the accretion rate onto the 
central black hole must be sufficient to power the luminosity 
of the AGN) . The characteristics of the disc should be de- 
termined by the rest of the nuclear system. Although /Edd 
changes on timescales much shorter than the evolution of 
the overall system, these are much longer than the dynamic 
timescales of the nISM. It is not currently feasible to deter- 
mine /Edd by modelling the disc and nISM simultaneously. 
We have chosen a simplified parameterisation so we can, to 
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first order, rescale /sdci, (see Section p.2[ ). Tlie rescaled solu- 
tions will have the same global properties as the real flows. 
Subsequent papers will fill in the important details in this 
picture. 

The flows in individual simulations are complex and 
varied, whilst exhibiting a generic bipolar circulation pat- 
tern. These results suggest why over 30 years of careful and 
detailed BELR observations and simulations have failed to 
identify a simple, global flow geometry. If anything, real 
flows will be more complex, not simpler, than the flows we 
find here. 

These generic bipolar flow patterns are characterised by 
both an accretion inflow and a wind outflow in the nISM. 
The gas which accretes can generate the luminosity assumed 
by the flow models. Bipolar outflows - and even hydrody- 
namical jets - are generated, so long as there is sufficient 
outward force along the accretion disc axis. These struc- 
tures will be found independent of whether the nuclei are 
radio-loud or radio-quiet: they are in general less coUimated 
and lower velocity than the most dramatic radio structures 
observed, although when a plume of gas escapes the nucleus 
it may be observable as a radio structure. 

The strength and structure of these flows depends pri- 
marily on the ratio of the gravitational and thermal energies 
in the nISM. Some of the inflowing gas is driven away into 
the wider galaxy, while some may accrete onto the black 
hole. The wind expelled from the nucleus can carry both 
significant mechanical energy and heavy elements into the 
wider galaxy. Thus, the nISM flows will affect significantly 
the structure of the ISM of the galaxy. At early times, the 
flows will enrich the galaxy ISM and may drive shocks and 



shells of gas out to distances of kiloparsecs (Beltrametti & 



1981 



Perryl j80|; Pyson, Falle fc Perry 198l| ; [Falle, Perry fc Dyson 



Ik fc Perry 1999) 



These nISM flows can only be observed indirectly. The 
nISM itself is hot - typically ~ 10^ K, the Compton temper- 
ature characteristic of the central radiation field - and does 
not give rise to observable line radiation, although limits on 
its column density are set by Thomson optical depths. The 
BELR is formed within the nISM, and so the broad line 
profiles will reflect the nISM velocities as well as the stel- 
lar orbital velocities. If the BELR is formed behind shocks 
driven by the stars as they move through the nISM (as in 
the models of PD) then the characteristic velocities of the 
line emitting material will be determined by stellar orbital, 
stellar mass-ejection and nISM velocities. The stellar veloc- 
ities will broaden the nISM contribution to the line widths, 
but the fundamental shape and asymmetries of the lines will 



reflect the nISM flows. E mpirical models of such profiles ( A.1 
bertsso i fc Perry 199S: ) show most of the characteristics of 
observed line profiles and their diversity. In this paper, we 
study the fiow structures which underly the structures de- 
termined by such empirical studies. 

We find that we can characterise the structure of the 
nuclear interstellar fiows on a two-dimensional plot. This 
is our first step in developing a theoretical diagnostic dia- 
gram for AGN flows, analogous to the Hertzsprung-Russell 
diagram which provides a systematic understanding of stel- 
lar structure. Stellar spectral type, luminosity, colour and 
temperature - characteristics which are determined by the 
basic physical parameters of age, metallicity and mass - are 
used as the axes of the HR diagram. Similarly, the axes of 



our plot can be expressed in a number of ways. Although it 
would be desirable to use the fundamental physical param- 
eters listed above as the axes, we find that the most natural 
variables are derived parameters which are closely related to 
the underlying physics (see Section 3.4). 



1.5 Scope and layout 

We present numerical models of cylindrically symmetric, 
isothermal fiows in black hole-starburst AGN. We have used 
global assumptions which are physically justifiable for QSOs. 
Various local deviations from those global approximations 
will be examined in detail in a future paper. 

In Section ^ we specify the parameters of the physical 
model of a supermassive black hole and a starburst stellar 
cluster. In Section ^, we give the hydrodynamical equations 
and discuss their scaling behaviour and their relationships 
to AGN observation. We also introduce the reduced set of 
dimensionless parameters with which we describe the gross 
features of the flows, to serve as a guide to the numeri- 
cal models which follow. The numerical methods used to 
implement the hydrodynamic representation of our physi- 
cal model are outlined in Section ^. We then describe the 
results of our detailed modelling. First, in Section |^ we de- 
scribe results for purely spherical flows. We verify that our 
methods reproduce the analytic results, and summarise the 
general properties of such flows, and those described in the 
literature. Next, in Section ^ we present the results from our 
aspherical models, calculated for parameters appropriate for 
luminous QSOs. In Section S.f: , we discuss models where the 



Eddington ratio is very small, as is suggested to occur in 
some Seyfert galaxies. In Section |^ we discuss these results 
and their physical implications in the context of the current 
observational and theoretical understanding of activity in 
galactic nuclei. We compare our results with the models of 
the flow structure in active nuclei of PD and Perry (1993a). 
We sketch the global constraints on the models. Readers 
who are most interested in the observational consequences 
of our work, and ways in which our model can be tested, 
could now turn straight to the flgures in Section ^ and then 
read Section ^, where we show that the symbiosis between 
a nuclear starburst stellar cluster and a supermassive ac- 
creting black hole, via the nISM, provides a self-consistent 
model for nuclear activity in galaxies. In Section ^, we draw 
together our conclusions. 

Illustrative animations are available at 



ht t p : / / ast . leeds . ac . uk/~ rj r w/agn . ht ml . 



2 PHYSICAL STRUCTURE OF THE ACTIVE 
REGION 

We model the nuclear system of an active galaxy as an ac- 
creting supermassive black hole, a dense starburst stellar 
cluster, and the resulting hot, flowing interstellar medium. 
In principle, mass, momentum and angular momentum will 
be exchanged across the interface between the active nu- 
clear region (radius ~ 1 pc) which we explicitly model, and 
the wider galaxy. In Section we discuss the theoretical 
and observational timescales which characterise our model. 
These determine the physical effects which dominate the 
structure of the hot phase flows. We define the parameters 
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Table 1. Glossary of non-standard symbols and abbreviations widely used in the text. 



Symbol Description 



BELR Broad emission line region 

Cs Sound speed in the nISM 

/(jisc Fraction of luminosity in accretion disc (anisotropic) component 

/Edd(^) Ratio of outward radiative forc e (f rom Thomson scattering) to inward gravity of the black hole 

{/Edd> Angle average of f-^dd, cf. eq. ( |a^ ) 

f(e) Form function for /Edd, so /Edd(S) = (/Edd>/(S) 

g, g Body-force acceleration (gravity and radiation) 

h Index of stellar distribution in stellar cluster halo, cf. cq. (^ 

H Thickness of accretion disc 

LboIi i'46 Accretion-generated bolometric luminosity of the BH, Lboi/10^^ ergs~^ 

L(jiac Accretion disc (anisotropic) component of LboI 

Lgdd Eddington luminosity of the black hole 

M Flow Mach number 

Ma, Rate of accretion of mass onto black hole 

Mc Mass of stellar cluster within rc (core mass) 

MciMclfi Total mass of the stellar cluster, Md/lO* Mq 

^h,Afh,8 Mass of the black hole, M© 

mi Low mass cut-off for stellar cluster IMF 

Mi, Mass of an individual star 

nISM Nuclear interstellar medium 

(r) Radial star number density of stellar cluster 

ric Star number density at edge of stellar cluster core 

n Total particle density 

Q, Q—s Rate of mass-loading from stellar cluster as a fraction of stellar density, Q/10~* yr~^ 

q{r) Mass loading rate 

(jc, q{r/rc) Mass loading rate in the core, variation of mass loading scaled to this value 

rc Radius of stellar cluster core 

Ta Sonic radius (for spherical models) 

iJ* Radius of an individual star 

R Radius of accretion disc 

s Index of stellar distribution in stellar cluster core 

tc Compton timescale at constant volume 

T4 Temperature T/10* K 

Tc,Tc,7 Compton temperature, Tc /Iff K 

Keplerian velocity (of circular stellar orbit) at cluster edge 

Vi, Escape speed from individual star 

r/acc Accretion efficiency factor 

A Dimensionless length scaling parameter 

/Xacc Fraction of stellar mass loss accreted by BH 

A** (f/rc ) Stellar ass enclosed with r as a fraction of total cluster mass 
p*(r), p*(r/rc) Radial mass density of stellar cluster 

p*_c Mass density of stellar cluster core 

a Velocity dispersion of the stars in the stellar cluster 

r Dimensionless time scaling parameter 

Tf nISM characteristic flow timescale, ~ lO-'^'^-lO^'^ s 

rcoii Stellar collision time 

Tdyu Cluster dynamical timescale, i.e. Tc/vk 

r* Characteristic stellar lifetime 

S Ionization parameter 

Dimensionless mass-input scaling parameter 



of our model in Section |2.2| , and describe our treatment of 
the interfaces between the flow and the accretion disc, a pos- 
sible nuclear jet and the surrounding galaxy in Section |2.3[ 
We provide a more detailed physical and observational back- 
ground to our assumptions in Appendix |^. 

Table |^ is a glossary of non-standard symbols and ab- 
breviations widely used in the text. 



2.1 Timescales 

We shall concentrate on global nuclear interstellar medium 
(nISM) flows, which have characteristic timescales Tf ~ 
10^°-10^^s, and on their response to global gravitational 
and radiative forces. In this regime, cooling timescales are 
short enough that the gas can be taken as in equilibrium. 
By comparison, relative timescales ^ 10* times larger are 
inferred for large-scale radio structures. In Table ^, we com- 
pare various physical timescales relevant to flows in AGN. 
Dynamical evolution timescales are long enough that 
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Table 2. Approximate physical timescales. 



Property 



Timescale (lO"' s) Scahng" 

Physical (A, r, 0) 



Stellar collisionf|^ 

Stellar relaxation in nucleus 

Galaxy interactions 

Stellar evolution 

Disc viscosity^ 

Orbital timescale 

Global sound crossing time 
Flow time 

Supernova evolution 
Post-shock coolinf^ 
Compton cooling 
Supernova frequency 
Cloud sound crossing 
Hot-phase ion collisions 
Cool-phase recombination 



3 X 10'^ 
1.6 X 10^ 
3 X 10^ 
3 X 10^5 
2 X 10-* 
30 

10 
1 

0.2 
0.2 
0.1 

0.03 
10-3 
10-4 
3 X IQ-^ 



Z/2 . ,-1/2 
Vc Afcl,8 

rpc(D/0.01c)-i T 
"r'^'(^SN/0.01c)-i T^/^r 



'pc^46 

Q~-lM-\ 



1/3 



A-2r3<j 
A-1t2 
A-3r2<; 



8-'"cl,8 
((i/10l3 cm)T4"^''^ 



° Variation of timescale with (A, r, 0), assuming all dimensionless parameters (e.g. nuclear spectrum, Eddington ratios) and 'microphys- 
ical' values (e.g. cool pha se tem perature) are constant. 

^ e.g. Perry & Williams (1993), assHHiing the stellar cluster dominates the gravitating mass. 

a ^ 1 is the Shakura-Sunvaev (1973) viscosity parameter, Td^4 is the temperature of the disc gas in 10^ K (see also Shlosman, Begelman 
fc Fran! 1990| ). 

Applies if post-shock cooling is dominated by bremsstrahlung cooling from 10^ K - for shock temperatures smaller than this, there is 
an extra factor T^/^ = A/r. 



the hydrodynamics is insensitive to secular changes in the 
gravitating mass distribution, although the details of the 
mas s loading do depend on the secu lar evolution of the clus- 
ter ( |MCD| purphy fc Perry 1999[ ). Stellar collisions, dy- 
namical relaxation, galaxy interactions, stellar evolution and 
disc viscosity all have characteristic timescales far longer 
than the flow time through the nucleus. These processes 
are all important in determining the secular evolution of 
the nucleus. Integrated over the nucleus, mass loss due to 
passive stellar evolution dominates that due to collisions 
for all but the densest (<: 10^ M0 pc"^) and most massive 
clusters. Mass loss due to relaxation of stars into the loss 
cone of plunging orbits are always insignificant. If the clus- 



gas will add to the disc at radii significantly outside 10"^ pc, 
so the feedback between fiow and central luminosity will 
have a very long delay (roughly 10^ times the flow time, 
Table This justifles our approach of calculating the flow 
structure for an applied central luminosity rather than at- 
tempting to find fully self-consistent solutions. 

Small timescale variability is most often measured in 
relatively low-luminosity Seyfert galaxies. Typical values are 
10~^Tf for continuum variability in X-rays, J; 10~'^Tf for the 
UV continuum and ^ lO^^rf for the optical continuum in 
Seyferts, which ar e both believed to b e driven by the high 
energy source (e.g. Edelson et al. 1996). Variatio ns of 10 per 



ter IM| extends well below livln, collisions can dominate 



cent in QSOs are found over timescales of 10 Tf ( Hook et al 



mass nipuL very cloHti to Lhti black hole (Murphy & Ptirry- 
199S|)~ Jence, if low luminosity AGN are characterised by 
small, centrally condensed starburst stellar clusters mixed in 



1994). In Seyfert galaxies, optical- UV emission lines also 
vary on a <: 10~*Tf timescale, with evidence for evolution of 
t he distribution of line-emitting gas over a lO'^Tf timescale 



dominate stellar 



a mass we old cluster, stellar collisions may 1 
evoluti( W aH the iiieclianiijui for tiiatjij Iohh iti riuch objec lij (cf. 



(^erry, van Groningen fc Wanders 1994 ; Wanders fc Peter-, 
son 1996|). The flow will be perturbed by sound waves driven 



MCD; WiUiams fc Perry 1994; Murphy fc Perry 1999) 



The similarity of the timescales for galaxy interactions 
and stellar evolution suggest that star formation will be im- 
portant through much of the life of the active nucleus. How- 
ever, threshold processes operating at intermediate scales 
within the galaxy may serve to reduce substantially the 
timescales which characterise the periods of mass input to 
the nucleus. The accretion disc viscous timescale is far longer 
than the dynamical timescale of the hot phase flow, even 
in the very central regions, ~ 10"^ pc, where temperatures 
reach 3 x 10* K (eqs A£ , AlC ) . We expect that the hot phase 



by short timescale variations in radiation forces. The steady 
luminosity we assume here may be taken as a time-average, 
which we assume drives a 'averaged' flow. 

The thermal evolution of the nISM may affect the flow 
at the 10 per cent level. This may be qualitatively impor- 
tant if instabilities result. The timescale for evolution of 
supernovae is comparable to the flow timescales, and so 
may be signiflcant - individual supernova events can sig- 
niflcant ly perturb the flow, particular ly in lower-luminosity 
nuclei (Perry, Williams & Dyson 199£). The short timescale 



for hot phase collisions is a necessary (if not sufficient) con- 
dition for our use of single-temperature hydrodynamics; the 
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timescales for sound to cross clouds and for recombination 
suggest that the hot phase may be taken as a steady back- 
ground in cloud structure problems. 

As we shall see (Section some flows seem to reach 
no equilibrium, undergoing a chaotic sequence of explosions 
or large magnitude oscillations; others seem to be steady 
on a dynamical timescale, but can take far longer than this 
to reach equilibrium. Flows with the latter behaviour can 
be modelled analytically in spherical symmetry (see Appen- 
dix ^). However, even those flows which do not reach equi- 
librium will adjust to changes in the structure of the nucleus, 
such as an increase in the central luminosity, within a flow 
crossing timescale. 



2.2 Components of Active Nuclei 

This paper concentrates on the structure of the hot global 
ISM which fills much of the volume of an active nucleus. 
In this section, we describe the other basic physical compo- 
nents of the nucleus which dominate its mass content and 
energetics, and which give rise to and act on the global flow 
to determine its structure. 



2.2.1 The Black Hole and Stellar Cluster 

We assume that an AGN contains a supermassive black hole, 
of mass Mh — lO*Mh,8M0, surrounded by a starburst stel- 
lar cluster, of total mass Md = 10*Afci,8 Mq. For the fidu- 
cial model of a luminous QSO . the cluster has a mass of 
Mci,8 ~ 10. Dynamical models (MCD) show that the stellar 



cluster dominates the gravitational potential of the nucleus, 
and that the black hole typically represents only 20 - 50 per 
cent of the mass of the syste m. O bservational support for 



this is discussed in appendices Al 



and A2 . Mass is liberated 



from the stars within the cluster and distributed through 
the volume of the nucleus. 

The stellar cluster is chosen to be steady and spherically 
symmetric with a stellar density given by 



■ p^.,cp*{r/rc), 



(1) 



where is the dimensionless spatial profile function nor- 
malised at the core radius, Vc. In this paper, we take a broken 
power- law form for (cf Appendix |A2|): 



(r/rc) r <rc 
(r/rc)"'' r > r-c 



(2) 



The power-law indices, s and h, are required to be s < 3 and 
ft > 3 for the cluster mass to be finite. In fact, we choose 
s = and ft = 5 in most of the models presented, similar 
to the Plummer model initial conditions of MCD. Exper- 
iments with different stellar distributions did not greatly 
alter the results that we discuss. Where we take ft = 5, the 
core radius is often chosen to be 2/3 pc, in order that the 
average radius of mass injection is 1 pc. With these param- 
eters, the central stellar density within Vc for our canonical 
model is p*,c = 3.2 x 10* Mq pc~^. This value is consider- 
ably larger tha n that inferre d from observations of nearby 
galactic nuclei ( Watson 199£ ) , but the lifetimes of the clus- 
ters we consider are short and their presence may be masked 
by the luminosity of the QSO they feed (for a more detailed 
discussion, see Appendix If the gravity and radiation 



forces are negligible in such a model, the outfiow becomes 
transonic at 0.89 pc. 



2.2.2 Stellar Mass Loss 

There are various sources of mass for the fiow in the nISM. 
Stars in the stellar cluster lose mass through the normal 
processes of stellar evolution, and also as a result of collisions 
within the stellar cluster and tidal disruption as they pass 
close to the black hole. The mass-loading rate, q, is chosen to 
be spherically symmetric, and time steady over the period 
of our calculations. Clearly the mass-loss rate evolves over 
the lifetime of the cluster and is a function of the age and 
IMF of the cluster. However, the timescales of this evolution 
are much longer than the typical dynamical timescales over 
which our calculations run. We write 



g(r) = qcq{r/rc). 



(3) 



In the present paper, we restrict ourselves to the case that 
the mass-loading rate to be proportional to the stellar den- 
sity, at a rate 



g(r) = Qp*(r), 

so that qc = Qp*,c. Q — 10~* yr~'^Q_8, and 
young cluster of massive stars (cf. A ppendix 
Perry 199^ [Murphy fc Perry 199£| ). 



(4) 
1 for a 



3; Williams 



2.2.3 Accretion Disc 

Accretion from the global flow, directly and through the 
accretion disc, powers emission from the vicinity of th e cen- 
tral black hole ( ^ycki, Collin- Souffrin fc Czerny 1995| ). The 
galaxy ISM may also feed an accretion disc if it does not 
hang up at an intermediate radius, and if t he extended 



Ipc scale) regions of the disc are stable (Hure et al 



1995) 



The accretion disc which forms around the black hole 
can also act as both a source and a sink of gas. In the hot 
central regions, a wind may be driven from the disc surface. 
Further out, gas may condense onto its surface, while self- 
gravitating instabilities could generate clumps of cool gas 
which would then be ablated by the nISM flow as it streams 
past. Eventually, the mass lost from the stellar cluster must 
either reach the black hole, to power the luminosity of the 
AGN, or be blown out of the nucleus in a wind. 



2.2.4 Luminosity 



The total integrated luminosity is 

r _ 1 n46 r -1 

LboI = 10 L46 erg s . 



(5) 



We model the angular distribution of the radiation field, 
Lsoiid), by two components, one isotropic - representing 
the radiation directly associated with the black hole - and 
the other anisotropic, representing the radiation of a thin 
accretion disc. This accretion disc luminosity is the only 
component of our model which is not spheri cally symmetric 



(Appendix A5). As we discuss in Appendix A6 



we assume 



for the moment that any variation in the continuum spec- 
trum with angle is unimportant. As we assume that the 
dominant opacity is electron scattering, we can describe the 
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Figure 1. Illustrations of the assumed structure of the nucleus. In a), a schematic view of the nucleus shows the extensive disc which 
surrounds the black hole. Along the axis of the disc, a radio jet is driven outwards in about 10 per cent of active nuclei. Around the 
black hole, the stellar density is high enough to be treated as a distinct stellar cluster - a small fraction of these stars are indicated by 
dots in this figure. The nebulosity which fills the region within the cluster illustrates the location of the nISM which we simulate. In b), 
we illustrate how these components relate to the structure of the simulations shown in the rest of this paper. The black hole is at the 
origin of the grid, the accretion disc in the horizontal plane, and any jet would travel along the vertical axis. The half-mass radius of the 
cluster reaches out a substantial way through the grid; mass loss from the cluster falls inwards in the plane of the accretion disc, and is 
driven outwards along its axis. 



luminosities by Eddington ratios (i.e., ratios of outward ra- 
diative force to the inward gravity of the black hole) for a 
given b [ack hole mass, where the radiation flux at (r, 9) is 



47rr^ 



(6) 



and where the Eddington limit luminosity (at which the ra- 
diation driving balances the gravitational acceleration) is 
given by 



LEdd ^ 1.3 X 10^^Mh,8ergs-\ 



(7) 



The angularly dependent Eddington ratio, /Edd(^), is spec- 
ified by its average over solid angle 



(./Edd) = ^ / ,fEdd(e)dn, 



(8) 



and an angular form function J{6) which depends on the 
fraction of the luminosity in the anisotro pic c omponent, /disc 
and its angular dependence (Appendix A5), so /Edd(^) = 
ifEdd) f{9)- Studies of AGN have indicated that typica l val 



ues of /Edd vary from 10 to values close to unity (Laor 



199C; C oUin-Souffrin 1993), with the higher ratios found for 
the brighter nuclei. However, even amongst Seyfert galaxies 
with very similar source luminosities the Ed dington ratio 
app ears to vary between lO"'^'^ (NGC 4258, [Lasota et al 



199e|)~alnd 0.5 (NGC 1068, [Greenhill et al. 199^ ) 



of the radiation, with values 0.6 ;$ ^ 15 predicted for 



the majori t y of AGN (Ikrolik, McKee fc Tarter 1981 



et al. 1986| [Mathews fc Ferland 1987] ) . We discuss the physi- 



Fabian 



cal limits on the validity of this assumption in Appendix A6 



The sound speed in fully ionized gas with solar abundance 
is Cs = 360r7^^ kms~^: we have used Cs = 10~"^r7'^^c for 
simplicity. 

The assumption of isothermality allows us to scale our 
models to a wide range of situations; however, there are some 
limits on the applicability of this assumption. The presence 
of strong shocks in our simulations means that including 
the energy equation explicitly would allow gas to cool: the 
production of such cool gas is important in deriving obser- 
vational diagnostics for the flow, in calculating the rate at 
which gas can be added to the accretion disc (though lim- 
its can be put on the amount of gas which can possibly be 
added to the disc in this manner), and in increasing the 
radiative force per unit mass, which may have important 
dynamical effects. However, in the context of global mod- 
els, the short time- and length-scales which will characterise 
this cool gas mean that it is difficult to accurately model 
both small and large-scale flows. Modelling the details of 
mass-injection from individual stellar sources presents sim- 
ilar difficulties. In consequence, we have decided to focus 
first on the global properties of flows which remain closely 
isothermal, and treat the finer-scale structures in future pa- 
pers. 



2.2.5 Thermal Equilibrium 

Consistent with our choice of electron scattering as the dom- 
inant opacity source, we assume that the gas is maintained 
everywhere in isothermal equilibrium, at the mean Comp- 
ton temperature of the central luminosity source, Tc = 
Tr X 10^ K. This temperature is determined by the spectrum 



2.2.6 Forces 

As well as hydrodynamic effects, body forces act on the 
nISM as a result of gravity, radiation pressure, and fric- 
tion with the mass injected into the flow. The gravity is 
dominated by the black hole and the stellar cluster (Appen- 
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dix . We assume that the mass loss from the cluster is, 
on average, at rest with respect to the black hole and stellar 
cluster and must therefore be accelerated by the flow. 

The radiative driving acceleration due to Thomson 
opacity is given by 



(9) 



where ctt is the electron scattering cross section. The grav- 
itational force on the nISM at r is due to both the mass of 
the black hole and the mass of the stellar cluster which lies 
within r; the latter can be written as 



Md(<r-) = AfciM*(r/rc), 



(10) 



where fi^{r/rc) is the fraction of the total cluster mass within 
r. T he f orm of /i* for our broken power-law is given in Sec- 
tion 



A4, 



The inverse-square laws for the central radiative and 
gravitational forces hold at all radii outside, respectively, 
~ 10~^ pc (where the finite spatial extent of the radiating 
accretion disc becomes important) and ~ 10~^ pc (where 
relativistic corrections need to be made), in QSOs. In ad- 
dition to these limits to the inverse-square law, the effects 
of the angular momentum of the hot phase, pressure from a 
thick central disc, and extra heating processes (such as stim- 
ulated Compton scattering) may effectively work counter to 
the gravitational force. We have chosen to mimic these ef- 
fects by smoothing the forces in a region close to the central 
black hole. Instead of a pure inverse-square law, we take the 
forces to vary as 



9 oc 



(r2 + e^f/^ ' 



(11) 



where in general we take the smoothing radius, e, to be 
far smaller than the cluster radius, but significantly larger 
than the grid resolution where possible. In this way, the grid 
resolves the fiow structures throughout the nucleus. 

We shall see that the link between the gravitational 
velocities at the radius at which the outflow is driven and 
the broad line kinematics leads to one of the fundamental 
distinctions between the results of the current paper and 
previous models of global flows in active nuclei. The low- 
ionization lines may be formed in exactly the regions of the 
accretion disc wh ich dominate the anisotropic continuum 
emission ([]DMP), which is one physical origin for such a 



smoothing term. This provides a link between LILs and the 
driving of the global wind in which the high-ionization lines 
are formed since these HIL are formed behind shocks within 
the global flow. The shock velocities, reflected in the line 
widths, represent the super position of the nISM velocities 
upon the stellar kinematics (Perry 1999). Both the low- and 
high-ionization emitting regions have velocities determined 
by the gravitational potential. The high-ionization emitting 
gas has an additional component of velocity due to the nISM 
flow. He nce the simil arity of high- and low-ionization line- 
widths (Corbin 1991) arises naturally in these models, al 



lowi ng for the centroid shift between the lines (Espey et al 
198S |; [l|ytler fc Fan 199^ ^ulentic et al. 1995[ ). 



2.3 Interfaces: Accretion discs, jets, the galaxy 

We now discuss the interfaces of the nuclear flow at the 
edges of our computational grid. In the plane z — 0, the 
flow interacts with the accretion disc. Along the grid axis, 
r = 0, there may be a relativistic jet. At large radii, the flow 
in the nucleus will interact with the galaxy ISM. 



2.3.1 Disc interface 

The nuclear starburst stellar cluster inevitably creates an 
nISM. No available treatment of accretion disc physics has 
considered the myriad physical processes which may arise 
in such an environment. The nISM flow near the disc plane 
can clearly have a wide range of detailed structures. At one 
extreme, the nISM and accretion flows are completely sep- 
arated by a contact discontinuity, and may be modelled by 
a mirror boundary condition. At the other extreme, if the 
pressure is sufficient that cooling can occur, all the material 
raining onto the disc will be lost to the global fiow and incor- 
porated into the disc. We model the latter case by specifying 
a region of the accretion disc over which we treat the disc 
plane as a mass sink. 

The scale height, H, oi & thin accretion disc is given by 



^ ^ 0.015 

r \ A/h.8 



1/2 



(12) 



where r4 ~ 1 is suggested by many models of the inner disc. 
This will be reduced still further by the gravity of the stellar 
cluster at radii beyond the cluster core. We can, at best, 
marginally resolve the cool disc in a (uniformly gridded) 
global model of the flow in AGN. The thermal equilibrium 
time of the disc gas is of order the recombination timescale, 
roughly 3 x 10* s, even at a density as low as 10* cm""^. To 
accurately model the structure of the disc would require a 
cell size ~ 10^* times the size of the global grid. 



2.3.2 Jet interface 

The other inner boundary condition on the models is along 
the axis of the disc. This will form the interface with a nu- 
clear jet, if one is present. While radio-loud AGN, with their 
strong jets, form only 10 per cent of the population, there is 
evidence for radio emission on BELR scales in many more 



nuclei (Dunlop et al. 1993; Kukula et al. 1997b) 



In the present paper, we treat the axis as a mirror 
boundary, in effect assuming there is no jet or that its width 
is ^ 1 pc. In reality, a boundary layer will form at the in- 
terface between global flow and jet. We leave the detailed 
treatment of such boundary layers for the present. In some 
cases we do find that well-coUimated non-relativistic jets 
form spontaneously in our model, driven by the wide open- 
ing angle outfiow from the central black hole. 



2.3.3 Galaxy interface 

In real active nuclei, mass, momentum and angular momen- 
tum will be exchanged between the active nuclear region 
and the wider galaxy. The rate of this exchange will depend 
upon both the properties of the nuclear region and of the 
galaxy. The timescales for stellar evolution and mass inflow 
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from galaxy interactions are comparable, but far longer than 
those over which the structure of the nuclear ISM is deter- 
mined. 

The rate of mass input by the galaxy to the nucleus has 
been the subject of much discussion. S hlosman et al. ( 1990, 
see also, e.g., the proceedings edited by Shlosman 1994) sug- 



gest that the rate may well be controlled by galactic inter- 
actions and subsequent instabilities. While some inflow of 
the hot phase from the ISM of a normal galaxy may oc- 
cur, it is difficult for this gas to lose its angular momentum 
and so the accretion rates will be a very small component 
of the mass budget of the nucleus (and may be further sup- 
pressed by any nuclear wind). Where a galaxy interaction 
has occurred, however, the angular momentum constraint is 
removed for gas distributed on kiloparsec scales. The simi- 
larity between interaction timescales and the stellar forma- 
tion and evolution timescales within the nucleus might then 
suggest that most active nuclei would be observed during 
the actively star forming phase. However, when the details 
of the fuelling process are considered, it is found that the 
infalling gas tends to remain well outside the nucleus un- 
til a sufficient mass has collected that it gas go through a 
self-gravitating instability and can then fall inwards once 
more (the so-called 'bars within bars' picture). The fuelling 
thus depends on this threshold process operating on scales 
of ~ 100 pc, where dynamical timescales are far shorter than 
those of galactic interactions. It is thus reasonable to con- 
centrate upon the case in which mass input from the wider 
galaxy is negligible in this first paper. 

Making only the reasonable assumption that the galaxy 
succeeds in creating a nuclear stellar cluster (cf. Watson 
1997, and references therein) we conclude that in situ feed- 
ing of the black hole by the stellar cluster is highly probable. 



3 THE THEORETICAL MODEL: TOWARDS 
AN HR DIAGRAM FOR AGN 

We next discuss the basic equations, moving on to the scal- 
ing of t he solutions to different sets of parameters (Sec- 
tion 3.2), the scalings which may characterise observed AGN 
(Section 3^), and the parameters which dominate the flow 
structures (Section 3.4, which also introduces our diagnostic 
plot). 

The interstellar medium is evolved using the equations 
of hydrodynamics, with a distributed mass-loading source 
function and gravitational and radiative forces. These are 
conservation of mass; 



| + V.(p.)=^, 

the equations of motion, 

Dv „ 

Pj^ = -Vp -qv- pg, 

and the energy equation. 



d (1 2 



= -V. 



, 1 2 

V \ -pv 



IP 



7-1 

-pg.v + p(r - nA) + 



(13) 



(14) 



(15) 



the mean mass per particle, 7 = 5/3 is the adiabatic con- 
stant, r and A are the usual heating and cooling coefficients, 
and Cq is the isothermal sound speed which characterises the 
internal energy of the mass injected. 

These equations include both the effects of applied grav- 
itational and radiative forces through the combined effective 
acceleration term. 



G[(l-/Edd(6l))Mh-HMci(<r)] 

(^2 ^ ^2)3/2 ' 



(16) 



where Mc\{<r) is given by equation (A6). Mass loading from 
the stellar cluster, at a rate q per unit volume per unit time, 
is included both in the mass conservation equation and as a 
frictional term in the equations of motion. The only radia- 
tive driving which we have included is that due to Thomson 
scattering. 



3.1 Isothermal Flows 

In the present paper we assume that the nISM is maintained 
in radiative equilibrium at the Compton temperature and 
hence do not solve the energy equation; instead we use the 
isothermal equation of state p — pc^ where the isothermal 
sound speed Cs is determined by Compton equilibrium (we 



discuss the validity of this assumption in Appendix A6) 



7-1' 

together with the equation of state p — pksT / p, where p is 



The input parameters to the equations governing the system 
are therefore the mass of the black hole Afh; the mass of 
the stellar cluster Afci, the radial distribution of the mean 
stellar density p^{r/rc), the cluster core radius r^; the mass 
loading rate density in the core q^, and its radial variation 
q{r/rc); the bolometric luminosity LboI, the sound speed in 
the nISM Cs and the distribution of the effective Eddington 
ratio /Edd(6'). 

Spherically symmetric, electron scattering driven, point 
source isothermal flows (/Edd(^) = const [less than unity]. 
Pi, = q = 0) have the well-known Parker solutions, with 
a characteristic length given by the sonic point radius, Vs, 
(rs = (1 — /Edd)Afh/Cs, Kippenhahn et al. 1975, hereafter 
KMP) and a characteristic velocity by the sound speed, Cs. 
In such flows, the radiation driving reduces the effective 
mass of the central gravitating object, Mh, to (1 — /Edd)Afh. 
As is well known, the solution (wind or accretion flow) which 
is realised in any particular physical situation is determined 
by the boundary conditions. 

This universal solution topology is removed by the ra- 
dial variation of the gravitational force or the distributed 
mass loading (both of which may introduce internal shocks 
and additional sonic point to the flow; mass loading also may 
introduce stagnation points), or, most dramatically, by the 
symmetry breaking introduced by the angular dependence 
of the radiation. Each set of /Edd(^), P*{'r'/rc) and q{r/rc) 
gives rise to a different family of solutions, further differ- 
entiated by the values of four dimensionless ratios between 
the other physical input parameters, as discussed below. Al- 
though this may, at first sight, appear to render any general 
conclusions impossible, the physical expectation that most 
AGN share similar spatial forms for /Edd, p* and q even 
when they differ in their total Mh, Md and LboI means that 
we can classify families of solution which apply over a broad 
range of objects, as we now discuss. In order to do so we 
write the equations in scaled, dimensionless form. 
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3.2 Dimensionless Equations 

We can cast equations ( p^ ) , ( p^ and (^6|) into dimensionless 
form by scaling all distances to the cluster core radius, rc 
(i.e. r — rc?, V — > (l/rc)V and e = rc?) and all velocities to 
the isothermal sound speed, Cs, (i.e. v{r,6,t) = Csv{r,8,t)), 
from which it follows that all times scale as t = rc/cgt. Fur- 
thermore p{r,9,t) = {qcrc/cs)p{r,8,t). We then have 

^^+V.(pi) = l (17) 

and 
Dv 

Pj^ = -Vp -qv- pg, (18) 
where 

_ GMci [(1 - (/Edd)/(g))(Mh/Mci) + M.(f)] 

^ rcCi (f2+g2)3/2 

Hence, from equations (^^, ( [l^ ) and (^^ we see that the 
dimensionless solution {v,p) depends on 

(l) f{0), p-i,{r) and q{f) - three dimensionless spatial form 
functions: these specify the angular variation of the radia- 
tion, the spatial variation of the stellar density and the mass 
loading, respectively, 

and the six physical parameters, Mh, Md, rc, LboI, Cs (or 
equivalently TnisM), and e, through the four dimensionless 
ratios: 

(ll) GMci/rcCg - essentially the ratio of the stellar cluster 
velocity dispersion to the sound speed, 

(ill) the ratio of the black hole to total cluster mass, Mh /Md , 
(iv) the Eddington ratio, (/Edd), and 
(v) the dimensionless smoothing length, e — e/rc. 

Note that no value dependent on q^ appears in the list, 
since the isothermal equations are homogeneous in p and q. 
Whereas qc sets the physical scale of p, it can be ignored for 
isothermal flows as far as the dynamics is concerned. This 
is a result of our assumption that /Edd(^) [or, equivalently, 
LsoiiS)] is a free parameter. If we had included the structure 
and dynamics of the accretion disc in the hydrodynamics 
then the rate of mass-loss, oc q, would determine I/boi(6') 
and /Edd (6') directly. 

By varying (l) through (v), each of the numerical solu- 
tions we present can be rescaled to a wide variety of AGN. In 
fact, as we shall discuss, two particular dimensionless ratios 
dominate the differences between solutions. The similarity of 
form of the solutions for a wide variety of AGN allows us to 
predict how, for example, the dynamics will scale with e.g. 
changing luminosity or Eddington ratio and spectral shape 
(or temperature) or luminosity and cluster core radius. 

For the solutions to scale, there are also constraints on 
the boundary conditions: for most of our simulations, these 
are observed as, for instance, free-flow boundary conditions 
which have no dimensional dependence; where accretion is 
allowed over some fraction of the disc, however, the size of 
this region must be specifled as a fraction of the cluster core 
radius for the models to scale. In reality, several of these 
dimensionless parameters will be less important than terms 
that break the scaling, such as heating and cooling of the 
nISM gas. 



As we shall see, within broadly similar bipolar nISM 
structures common to all symbiotic AGN, there are system- 
atic differences in topology (see Fig. ^). AGN of widely dif- 
ferent bolometric luminosities can have topologically iden- 
tical flows if they all share a common stellar density profile 
(i.e. core halo structure p(r)), have stellar mass loss rates 
which are directly proportional to p*, have the same radia- 
tion pattern f{6), ratios (ll) through (v), and scaled bound- 
ary conditions. I.e. if their underlying physical structure is 
similar, their nISM will be described by a single dimension- 
less solution (i3, p) to equations (|l^), dig ) and ([l9|). However, 
rather than present dimensionless solutions for particular 
combinations of (l) through (v), we have chosen to present 
all our results in physical units for fiducial models in order to 
make them more immediately accessible. Below, we describe 
the manner in which the models which we have calculated 
can be scaled to a wide range of parameter values, and dis- 
cuss how these scalings may relate to the relationships which 
have been inferred observationally for AGN. 

We have grouped p*(f) and q(f) in (l) above, since in 
this paper we make the further restriction that these spatial 
distribution functions are equal, which is equivalent to the 
assumption that qc = Qp*,c- [For this limit to be strictly 
valid, we must consider young stellar clusters with core den- 
sities pt,c ^ 10® Mq pc~^, above which collisions and tidal 
disruption at s mall ra dii make a significant contribution to 
the mass loss (MCD). Models A, B, C are all below this 
limit, but may breach it when scaled to Seyfert parameters.] 
In this case equation dlTI) becomes 



dp 



dt 



= -f- V. (pv) 



(20) 



The densities and pressures in the nISM flow are then p — 
{Qp*,crc/cs)p and p = (Qp*,crcCs)p. 

Clearly in such isothermal flows, a consequence of the 
scaling is that changes in the source spectrum, which result 
in a change in the Compton temperature and thus in Cs, 
imply a consequent change in the time scales. Other things 
being equal, the same t corresponds to a shorter real time 
in hotter flows (with harder spectra). Note, again, that no 
value dependent on Q appears in the list. 

Taken together, if ratios (ll), (ill), and (iv) are constant, 
then each of LboI, Afh and Md is proportional to rcCg with 
a known factor, so LboI oc Mh oc Md. 



3.3 From High Luminosities to Low: How AGN 
Scale in Nature 

One of the most striking features observed in the spectra 
of AGN is the remarkable consistency of their dynamics 
over orders of magnitude in total bolometric luminosity. Line 
widths - direct measures of gas velocities - vary almost im- 
perceptibly as Z/Boi changes from < 10'''^ to > 10'*'^ ergs" ^. 
The size of the emitting region, however, appears to track 
the luminosity, rBLR oc L^^^ (from reverberation and pho- 
toionization studies). Controversy still surrounds the ques- 
tion of whether or not black hole masses and/or the effi- 
ciency of the radiation process ((/Edd)) vary as LboI varies 
- or, more accurately put, whether LboI varies because of 
variations in Mh, /Edd, or because of the availability of fuel. 
The fuel supply is controlled by stellar and galactic pro- 
cesses. Unfortunately, to date observations can measure nei- 
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ther the properties (Afd, p{r), age or IMF) of the stellar 
cluster (because it i s outshone by the cent ral radiation by 
orders of magnitude Williams & Perry 1994), nor the rate of 
any possible inflow of matter from the galaxy. The relation- 
ship between the properties of the cluster and the fuelling 
process is, at least in part, a theoretical prediction of this 
paper. 

We now consider the implications for AGN of the scale 
invariance of the equations and discuss how the physical 
properties of AGN are expected to vary as a function of 
luminosity, I/boI- Our models make clear predictions about 
how unobserved properties (e.g. Mc\ and r^) are expected to 
behave as a consequence of those which are observed. They 
also provide simple rules for scaling the solutions we present 
here to objects with other luminosities or cluster/BH prop- 
erties so that the reader can choose amongst the solutions we 
present those which are appropriate to her or his favourite 
AGN. 

Consider a family of structurally similar AGN (i.e. the 
same (l)) all of which share the same characteristic values 
of (ll), (ill) and (v), but which differ in LboI, and possibly 
in (/Edd) (iv). The members of the family will have differ- 
ent characteristic sizes, velocities, time scales and masses. If 
their spectral hardness differs, so will their Compton tem- 
peratures, and thus their sound speeds. We denote the ra- 
tio of their sound speeds by v, the ratio of their typical 
length scales by A, and the ratio of their timescales by 
r. For consistency, r oc X/i/. Further, we denote the ra- 
tio of their mass loss rates by (j) (which, for given A and 
ror v, sets the gas density scale). The constancy of (ll) 
means that cluster masses scale as Md cx the stellar 

mass densities as p*,c oc {y/Xf and the nISM gas den- 
sities as p oc <j)u/\ oc (ji/r. The ionization parameter in 
the flow, H — Li/ (A-Kr^nnkTc) (see Section 7. J), scales as 
H oc L/4>\V^. 

3.3.1 AGN with constant Eddington ratios 

We first consider a class of AGN of different total lumi- 
nosities which are all able to accrete with equal efficiencies 
and, in consequence, to radiate with the same constant Ed- 
dington ratio. In such a population, with (l) through (v) 
constant, the masses of the black hole and stellar cluster 
scale as the luminosity (Md oc Mh oc I//(/Edd) oc L) and 
oc A^/r^ oc L; the ionization parameter then scales as 
S oc l/(j)u. 

(a) If, in such an AGN population, the velocities were con- 
stant (roughly as observed), the characteristic sizes of the 
regions would scale as A oc L (rather than A^ oc L as is often 
discussed), and typical timescales would scale as r oc A oc L. 
The stellar clusters would be very much denser and smaller 
as the luminosity decreased: A oc L, p*,c oc Further- 
more, the ionization parameter would be a function only of 
the mass loss rate of the cluster stellar populations through 
(/): E oc l/4>. 

(b) Alternatively, \i L/r^ were constant within a p opulation, 
as is often inferred from broad line studies, (e.g. Netzer & 



Peterson 1997), then A oc L oc Mh. In that case, our models 



1/2 



predict that the velocities would scale as oc (L/X) 
L^/"*, implying also that the spectra of AGN must harden 
slightly {{hu) oc T oc L^^^) as they become brighter. The 
stellar clusters would be more compact and denser as the 



luminosity decreased, but the changes would be significantly 
weaker than in (a), A oc L^^^, p^^c oc which seems 

more probable. Furthermore, the typical timescales and the 
ionization parameter would be weakly dependent on L, r oc 
L^/" and Hoc 1/0L^/^ 

There is some observational evidence for this kind of re- 
lation betwe en lin e width and luminosity. Baldwin, Wampler 
& Gaskell (1989) find a correlation between the width of 
C III] 1909 and source luminosity, but not for the other 
lines they studied. The variation in Compton tempera- 
ture also compares well with obser ved change in continuum 
shape with source luminosity (e.g. Zheng & Malkan 1993). 
Mushotzky & Ferland (1984) find an anticorrelation between 
ionization parameter and luminosity of this form can also 
explain the Baldwin effect, an anticorrelation between the 
equivalent width of lines and the source luminosity. All such 
observations are, however, plagued by selection effects and 
the wide dispersion between the emission line properties of 
individual galaxies with a certain luminosity. 



3.3.2 Eddington ratio varying with LboI- 

It is possible that the Eddington ratio varies systematically 
with source luminosity (a point to which we return below). 
This would seem to exclude scaling the solution to equa- 
tions (jl|) and ^) found for one AGN to other AGNs with 
different (/Edd) since, as we have just seen, (/Edd) is one of 
the dimensionless ratios defining the topology of the solution 
which we wish to apply to AGN of varying LboI . However, as 
we shall see from our simulations, the most important value 
of the Eddington ratio for the flow is that on the axis of the 
accretion disc, /Edd,o- Thus, for the moment, assume that we 
can ignore the angular variation of the luminosity. In that 
case, we can relax the condition that f{6) must be invariant 
(this assumption is supported by our numerical simulations, 
see Section o) , which is equivalent to allowing the degree of 
beaming to vary with L or /Edd- Since theoretical studies 
of the structure of accretion discs find structural changes 
with /Edd, this seems a reasonable assumption. When f{9) 
is constant (/Edd) oc /Edd.o, whereas when f{8) varies with 
L, (/Edd) and /Edd,o are independent. 

Given the primacy of /Edd,o, we find that the forces 
which act upon the flow are controlled by an effective black 
hole mass, (1 — /Edd,o)Afh, rather than directly by f{0) and 
the BH mass separately. Parameters (ill) and (iv) then com- 
bine into a single parameter, (1 — /Edd,o)(-^^h/A^ci) which 
defines the solution and must be scaled for different AGN. 
Note that in the trivial case that the radiation is isotropic 
and /Edd(^) = 1, the BH has no net effect on the flow and 
its mass is arbitrary. 

Consider now such a class of AGN whose Eddington ra- 
tios, (/Edd), vary systematically with luminosity. Requiring 
(./Edd.o — l)Mh/Mci to be constant rather than (ill) and (iv) 
introduces an additional degree of freedom, and relaxes the 
strict scaling between the masses and the luminosity. The 
constancy of (ll) means that, as before, Md oc . Since 
Mh oc L/(/Edd), the new condition, (/Edd.o - l)A/h/Md = 
const., becomes (/Edd.o - l)/(/Edd) oc Xu^/L. 
(c) Consider first a population for which the velocities vary 
as weakly as i/ oc L^^^, for which (/Edd.o — l)/{/Edd) oc 
X/L^'''^. If, in addition, L oc A^, (/Edd,o — l)/{/Edd) = const.; 
as we found in (b) above, v oc L^'''* combined with L on 
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implies that /Edd is constant; all other properties then vary 
as they did in (b). 

(d) If the velocities have constant characteristic values [v = 
constant), then the masses of the clusters must scale as 
Mc\ oc A whereas the masses of the black holes scale as 
Mh oc L/{fEdd} oc A/(/Edd,o - !)• If, in addition, L oc A^, 
we find that the Eddington ratios must scale as (/Edd.o — 
l)/(/Edd) oc L~^^^ . Typical timescales, r oc A oc L^^^, in ac- 
cord with the observation that typical timescales are shorter 
in low luminosity objects. We find also that H oc L^^^/(j>, 
i.e. for clusters with the same IMF and age, {(j) ~ const), 
HocLi/^ 

If the relation (/Edd.o — l)/{/Edd) oc L~^^^ found in (d) 
is imposed, and we attempt to impose f{6) constant, we find 
that both /Edd.o and (/Edd) increase as L decreases. This 
correlation is in the opposite direction from that inferred 
from observation, but will occur only for a narrow range 
of luminosity before the luminosity becomes super Edding- 
ton at all angles. The flow structures change dramatically 
when this occurs, suggesting that they will not scale be- 
tween Seyfert galaxies and QSOs. Alternatively, the flows 
can scale when the Eddington ratio varies with L and the 
characteristic velocities are constant, but only if the contin- 
uum becomes closely beamed as the luminosity decreases. 
If the on-axis Eddington ratio is roughly constant, we find 
that (/Edd) oc L^^^ which accords well with some estimates 
of the variation of the variation of /Edd with L. Furthermore, 
this then means that Mh oc A/ci oc L^^^; the stellar cluster 
will be more compact and denser as L decreases, A oc L^^^ , 
p^,c oc 1/L. This situation seems in broad accord with ob- 
servations. We conclude that in the case of constant /Edd.o 
but varying (/Edd) we can scale our solutions for QSOs to 
low-L AGN. 

The variation of some physical timescales with these 
scaling parameters is included in Table ^. If dynamical ve- 
locities are constant, most of the timescales will vary in pro- 
portion to the length scale; QSO models can then be scaled 
to other situations. Notable exceptions are the stellar col- 
lision time, and the frequency and lifetime of supernovae. 
In smaller nuclei, stellar collisions will become an important 
mass source while supernovae will become infrequent events: 
when they do occur, both processes will have a catastrophic 
effect on the nuclear ISM. 



3.4 An HR Diagram for AGN? 

As we have discussed, our models are determined by four 
dimensionless ratios and a number of 'form factors'. How- 
ever, extensive numerical modelling has suggested that the 
number of parameters which dominate the overall form of 
the flow is far smaller than is yielded by such a formal re- 
duction. We find that the solutions we calculate can be well 
described by just two dimensionless ratios, describing how 
strongly the gas is bound to the nucleus as a whole and how 
strongly it is driven away by the black hole at its centre. In 
the rest of this paper, we present models to illustrate the 
reasonableness of this simplified parameterisation, and its 
usefulness as a classification scheme for the numerical mod- 
els. The details of this classification scheme are presented 
in Sections 3.4.2 and While the derivation of detailed 



Table 3. Derived velocities, for the models defined in Table Q 
The values of vj^, Vz. and Cs are set by the parameters of the 
model. The ejection velocity, Vf,j , is the peak outflow speed in a 
region close to 2; = 0.2pc on axis (except f E close to z = 0.04 pc, 
F close to z = 2pc), where the mass ejected from the region close 
to the black hole has reached close to its final value. 
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Figure 3. In this plo t, we compare the outflow velocity, es- 
timated from eq. (Elh with that measured in the simulations, 
Voj ■ Values are given in Table ^. A good relationship is seen, 
apart frorn in Models E and F where the assumptions under- 
lying eq. ( |2l| ) break down. In Model E the smoothing region flUs 
much of the volume of the cluster so the cluster mass acts to sup- 
press the outflow. In Model F, the smoothing region is unresolved 
by the grid, so the efi'ective e will be larger than that imposed. 
Models B, D and H are coincident. 



observational diagnostics is beyond the scope of the present 
paper, we intend to work to bridge the gap between em- 
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Clegg (1985; cf. also KviUiams & Dyson 19941); 



V ./ c 

Figure 2. This graph compares the values of tik/cs and v^j/cs, also presented in Table 
to various models; CC - the uniformly mass-loa ded, z ero-gravity models of Chevalier 
BMS — the disc wind models of Begelman et al. (1983); KMP - the nuclear wind models of Kippenhahn et al. and Beltrametti & Perry 
approximate the on-axis flows here. The small panels show the outline structures expected for each set of parameters: arrows show the 
directions of supersonic flows, solid lines show the most important shocks and dashed lines the most significant sub- to super-sonic 
transitions. The regions of parameter space in which we find explosions and coUimated jets are labelled. The set of arrows to the right 
shows the directions in which models change if the variable labelled is changed. 



pirical classifications and this theoretical framework. Before 
describing the classification scheme we preview the physical 
properties of the flows we find. 



3.4-1 A preview of the nISM flows 



Here, we preview our results, and move on, in Section 3.4.2 
to the classiflcation scheme which we have found. The details 
of the models are presented in Section By anticipating the 
conclusions here, we hope the reader will be better placed 
to understand not only the details of the individual models, 
but also their relative place in the wider scheme. 

In spherical models with no accreting central black hole 
the form of the flow is principally determined by the ratio 
vk/cs of the Kepler velocity and the sound speed in the ISM 
(see also Durisen & Burns 1981), where «k = GMc/tc and 
Mc is the core mass of the cluster. For small vk ^ Ca, the 
cluster core fills with gas on a dynamical timescale and the 
sonic surface coincides with the edge of the cluster. By con- 
trast, for large vk/cb (as in Fig. most of the gas injected 
falls supersonically through the cluster onto a hydrostatic 
core at early times. Eventually this hydrostatic core fills up 
to the edge of the cluster, causing the infiow to cease. After 
this the core grows further, quasi-statically, until it reaches 
the sonic surface at rs ~ (i'K/cs)'^rc, when a steady isother- 
mal wind is established. This scheme is illustrated by the 
analysis of Appendix ^ and the numerical models in Sec- 
tion |5[ 

The addition of an accreting black hole to the cluster 
fundamentally alters the structure of the nISM created by 



the stellar mass-loss. The non-spherical radiation, by exert- 
ing a non-spherical force on the gas, gives rise to the in- 
and out-flow structures which we discuss. We have consid- 
ered a wide range of models spanning the expected range of 
the independent parameters listed above. In this paper we 
present a selected set of characteristic models which illus- 
trate the full variety of expected flow topologies. Our basic 
model, designated A, is the most likely physical scenario. 
We consider the systematic effects of changing the physi- 
cal structure of the region: in models Band C we decrease 
the relative mass of the cluster; models A^c explore accretion 
(where x indicates the region of accretion); D and E have 
a dynamically less significant cluster; F has a diffuse cluster 
and a heavy black hole; G has a small and diffuse cluster. 
The properties of our models are given in Table ^. 

In these two-dimensional models, gas which falls close 
to the black hole is driven outwards along the axis in a dy- 
namical timescale. Rather than a hydrostatic core forming 
in the gravitational well of the central black hole, a circula- 
tory fiow forms close to the centre. In addition, the amount 
of mass within the nuclei can continue to vary, either be- 
cause no gas escapes, as in Model C, or because it escapes 
in an episodic fashion, as in Model B. 

The properties of the flow may also depend on the age 
of the nucleus. Where the central outflow drives directly out 
into the ambient ISM, the nuclear flow rapidly reaches an 
equilibrium situation: it is unlikely that the transient evolu- 
tion, over a dynamical timescale, will be observed (or even 
occur, depending as it does on the artificial initial condi- 
tions). Age will certainly be important where most of the 
mass input is retained by the nucleus, since the mass con- 
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tent of the ISM will then increase linearly with time. Even- 
tually, this will lead the flow to break our assumptions of 
low optical depth and isothermality. In the mean time, as 
the mass content increases, the dynamical importance of the 
mass loading decreases, so the density scaling works for us 
once more. The structure of the flows becomes a turbulent- 
looking plume, with an extent determined by the balance 
between driving and gravitational velocities, in which the 
gas density increases linearly with time, surrounded by a 
constant accretion flow which is no longer dynamically im- 
portant (or observationally, as a result of its low gas pres- 
sure) . 



3.4-2 The classification scheme: An HR diagram for 
AGN? 

Many authors have attempted to find out how many ef- 
fective parameters determine the structure of active nuclei, 
using both observationa l data and theoretica l analysis (e.g 



Boroson & Green 1992 



Yi 199e 



Francis et al. 1992 



Blackman & 



Donea & Biermann 1996). This may be seen as a 



search for diagnostics analogous to the Hertzsprung-Russell 
diagram, which has led to so many insights into the most 
stars very close to spherical symmetry: it is no surprise that 
the situation in AGN is substantially more difficult. 

To recap, we take as the starting point for our model a 
galactic nucleus with an accreting black hole, surrounded by 
a starburst stellar cluster. We assume that the cluster can 
form by prior processes such as galaxy interactions. These 
two simple, relatively well-understood input components 
turn out to provide a remarkable and fascinating variety 
in the behaviour of the nISM. Broadly speaking, the results 
of our simulations can be understood in a two-dimensional 
theoretical classification scheme, analogous to the HR dia- 
gram. 

We find that we can classify the models computed in 
terms of two dimensionless quantities: vk/cs, the ratio be- 
tween stellar dynamical velocity and ISM sound speed; and 
Vcj/cs, the Mach number of the gas driven away from the 
black hole along the axis of the disc. These parameters are 
the axes of in Fig. ^ Also shown on this figure are the codes 
of the various models which will be presented in Section ^ 
and the form of the fiows which characterise each region 
of the plot. The flows take a variety of forms: inflow with a 
small core, jets, explosions and bipolar winds, as will be seen 
from the models themselves. The application of this plo t to 



classifying AGN models is discussed in detail in Section 7.1 



The interaction of the axial outflow with the gravity 
fleld and the mass input of the stellar cluster determines 
the form of the global flow: the axial wind can be confined 
within the core, can escape episodically or can drive out- 
wards as a steady flow. The overall effect is determined by 
a combination of whether the central outflow has more than 
escape velocity from the nucleus as a whole (with due ac- 
count taken for mass loading), and, if so, whether it has 
sufficient momentum to overcome the ram pressure of the 
accretion fiow. 

The most important effect not included in the models 
presented here is the thermal balance of the flow. As we as- 
sume that the flow is isothermal, the two-dimensional clas- 
sification of disc winds introduc ed bv Begelman et al. (1983; 
and related to global flows by Williams 1999 ) collapses to 



the axis labelled BMS on Fig. ^ - the one-dimensional clas- 
sification of Section ^ While isothermality seems to be a 
reasonable first approximation for flows in AGN (Appen- 
dix A6), we intend to investigate the structure of radia- 



tive flows in more detail in future papers. This classiflcation 
must, of course, be supplemented by the inclination angle to 
the observer when applied to observations. 



4 NUMERICAL METHODS 

These simulations were compu ted with a Godunov-type hy- 
drodynamic code (Falle 1991). The method is second or- 
der in space and time in smooth regions, representing the 
flow with piecewise-linear variations across the computa- 
tional cells. It calculates the inter-cell fluxes by solving a 
Riemann shock-tube problem across the cell interfaces. In 
the models presented here, the flow is treated as isothermal. 
The computational grid is two dimensional, with interfaces 
along the ordinates of a cylindrical polar grid. 

we assume that no mass 



As discussed in Section 2.3 



passes inwards from the galaxy to the nucleus at the outer 
boundaries of the grid, at large r and large z. If the nuclear 
flow at the grid boundary is outwards we use zero transverse 
gradient boundary conditions, and keep a record of the total 
mass and energy flows through these boundaries. This allows 
us both to verify mass conservation and to investigate the 
strength of the winds interacting with the external medium. 

In some models, the mass escaping the region has too 
little kinetic energy to escape to infinity. It is likely that some 
fraction of this mass ought eventually fall back into the nu- 
cleus, probably in the accretion disc plane, after some inter- 
val. We comment elsewhere (Appendix A2) on the insensi- 
tivity of our results to the detailed form of the mass-loading 
(whether internal or extra-nuclear) if this input mass has low 
angular momentum. In the models presented here, we have 
chosen to keep an independent total of the 'bound ejecta' to 
compare with other mass fiuxes rather than impose an arbi- 
trary reinjection model. Experiments with larger grids and 
a variety of reinjection models have shown that while this 
artificial lowering of the barrier to escape from the nucleus 
changes the boundaries between the flow modes by a small 
amount, it does not fundamentally alter the nature of the 
flows. 

We do not include the effects of angular velocities about 
the grid axis. We consider only the input of gas with low net 
angular momentum from stellar mass loss. Models including 
disc winds, although not the gravitational effec t of a stellar 
cluster, have been calculated by Woods et al. (1996). They 



require treatment of both angular momentum and of wind 
acceleration very close to the accretion disc. The z com- 
ponent of angular momentum is included in those models. 
While essential in that case as the mass source had high an- 
gular momentum, this results in nearly empty vortices form- 
ing along the poles of the simulation, which may not be phys- 
ically realistic when account is taken of non-axisymmetric 
perturbations to the fiow. 

We concentrate here on two limits, one in which there 
is negligible net mass transfer to the accretion disc (so that 
the surface z = is just treated as a mirror plane), and 
another in which mass is accreted over some fraction of its 
surface. We determine the flux onto the accretion disc in the 
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latter case by making the disc surface a free flow boundary. 
This mimics the effect of the strong cooling which may well 
be driven by the high pressures found in the disc plane (e.g. 
Figs. I @). 



5 SPHERICAL MODELS 

One limiting case of the flow models presented in this paper 
is the flow in a purely spherical system, such as an active 
nucleus with no anisotropic radiation, or a starburst with 
no black hole, or a globular cluster. Prior to presenting our 
main results, we present one example of a spherical flow as 
a test case. This can be compared with many treatments in 
the literature, with a variety of distributions of gravitating 



Table 4. Computational Grids 



matt e r and radiative forces (e.g. KMP ; Beltrametti & Perry 
198C|; iDavid, Durisen fc Cohn 1987] ). We focus on a case 



with no black hole and no radiation forces to illustrate the 
accuracy of the numerical method. Our numerical results 
agree well with an analytic treatment of the flow structure. 
In Fig. |4|, we show the structure of the flow after 



4.4 X 1 ) yr. An accretion shock (0.5 pc from the centre of 
the grid) then surrounds a central hydrostatic core, and is 
still moving outwards. No physical instabilities operate to 
drive the flow far from spherical symmetry. The structure is 
that of a central, near-hydrostatic core surrounded by a su- 
personic accretion flow. The velocities in the core (Fig. ^) 
reach levels as high as Mach 0.7 in one isolated region, and 
some striated structures appear close to the bounding shock. 
Some turbulence in the core would be expected in a real 
cluster, as a result of stirring by orbiting stars and of inho- 
mogeneities in the infalling flow. Here, however, turbulence 
is driven by the jagged form of the well-resolved shock at the 
cluster edge and by the well-kno wn numerical instability of 
slow- moving shocks ( Cjuirk 1994 ) , and is maintained in part 
by the reflective boundaries at r = and z = (notice the 
vortex at 0.15, 0.6; cf. the similar effects seen by Mellema & 
Frank 1995, for outgoing flows). 

Beyond the sharp edge of the cluster, a near-hydrostatic 
atmosphere reaches outwards. The flow would only become 
transonic at a radius of about 25 pc (i.e. 2.5M8r:^^ pc), far 
beyond our grid. The steep density gradient in the atmo- 
sphere allows us to ignore these outer regions as they will 
have little dynamic effect, and only a negligible influence 
on the overall mass budget. This density gradient amplifles 
small waves at the edge of the core to generate sound waves 
with gas velocities up to Mach 0.2 further out in the halo. 
Nowhere, however, are the density contours significantly per- 
turbed from spherical symmetry. 

In Appendix ^, we give an analytic derivation of the 
time evolution of the mass-loaded isothermal flow in an ini- 
tially evacuated spherical cluster with no central black hole. 
In Fig. 



Bl 



we compare these results to the numerical model 
shown here, and flnd excellent agreement. The smooth curve 
on Fig. Bla shows the good fit of the density to the isother- 



mal, hydrostatic solution obtained below (equation B9). The 



low level features of Fig. Bl d within the hydrostatic core are 
caused by the turbulence in the core. 

As the flow evolves, the hydrostatic core will become 
larger until eventually its outer edge, having merged with 
the inner accretion sonic point, reaches the outer sonic point 
of the flow and forms a classic Parker-type wind. Even after 



Models Colls Sizo/pc Smoothing radius/pc 



0ABCA,D 


300 


X 


300 


1.5 


X 


1.5 


0.05 


ESaSn 


300 


X 


300 


0.15 


X 


0.15 


0.05 


F 


300 


X 


300 


15 


X 


15 


0.05 


G 


1000 


X 


1000 


4 


X 


4 


0.004 


H 


600 


X 


600 


1.5 


X 


1.5 


0.05 



more than 4 x 10 yr the flow is far from equilibrium (cf. 
Figs, y and Bl). It is likely that by the time the equilib- 



rium solution is reached the stellar cluster will have evolved 
signiflcantly. As we discuss below, the high densities which 
are reached in the hydrostatic core are also likely to allow 
bremsstrahlung cooling to become important before the core 
has flUed, breaking our assumption of isothermality. 

Adding a black hole with no aspherical radiation com- 
ponent or accretion simply changes the properties of the 
central hydrostatic core. Small accretion r ates can be mod- 
elled as Bondi accretion within this core (Durisen & Burns 



1981; David, Durisen & Cohn 1987): as the accretion sonic 
point moves further out into the core, the mass removed 
will weaken the wind which eventually drives outwards from 
the cluster (cf. Model Sa below). In the next Section, we 
illustrate the dramatic effects which can occur when the ra- 
diation field of the central black hole and accretion structure 
is aspherical. 



6 ASPHERICAL MODELS 

We present in this section the results for a variety of models 
chosen to illustrate the relative importance of the black hole 
mass and luminosity, the nuclear cluster and the structure 
of the accretion disc. If the outward driving by the accretion 
disc luminosity can counteract the gravitational acceleration 
of the black hole, its anisotropy inevitably results in a bipo- 
lar outflow. This outflow coexists with inflow from the outer 
regions of the cluster. The outflow funnels the inflow close 
to the accretion disc plane at small radii. 

We explore a range of models between those in which 
the central flow is strong enough to drive a steady outflow, 
with a classic 'bipolar' shape, and those where it is confined 
(smothered) within a very small central region. At inter- 
mediate values, the fiow either drives jets or sequences of 
explosions, depending on the ratio of the Keplerian veloc- 
ity to hot phase sound speed. The balance between infiow 
and outflow is determined primarily by the ratio between 
the dynamical velocity in the cluster core and that at the 
smoothing radius around the central black hole. 

The generalizations we draw throughout this paper 
about the pattern of structures which occur are based on 
numerous simulations - of these, those which we have chosen 
to present here are those with parameters close to those ex- 
pected from the observations (as discussed in Appendix^, 
and those where the resulting nISM structures are non- 
trivial. If, for instance, the stellar cluster has a mass far 
lower than the central black hole, then its influence on the 
flow will be small, a situation already treated thoroughly in 
the literature. 

In Tables ^ and ^ we list the parameters of the mod- 
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a) b) 

Danaity/em"* TalMity/O.Olt: 




r/pc r/pe 

Figure 4. A spherical flow model (Model discussed in in Section ^, with parameters given in Table |5|) at 4 = 4.4 X 10* yr. The panels 
show a) number density n (up to 5.6 X 10^ cm~^) b) total velocity u (grey-scale and contours, up to 3.7 X 10~^c). The vectors are in the 
flow direction, and have lengths proportional to the velocity, sampled at one vector in every 15 X 15 cells. The flow is sonic at the w = 0.1 
contour, c) and d) are a magnification of the central 0.5 X 0.5 pc region (here, the vectors are every 5x5 cells). The flow structure is of 
spherical accretion onto a near-hydrostatic core (Appendix nBl). 



els presented, and in Table ^ various numerical results. As 
can be seen, we first explore (in Models A, B, C) the bal- 
ance between the dynamical velocities in the cluster and at 
the smoothing radius by varying the black hole parameters, 
leaving the cluster parameters unchanged. Next, we present 
models in which the accretion disc is assumed to remove gas 
from the flow over some fraction of its surface (Ao.oos, A0.25, 
-^'0.25)- We then investigate the efTect of varying the param- 
eters of the stellar cluster in Models D, E and F. We look 
at the effects of changing other flow parameters, such as the 
opening angle of the central outflow (H), and of increasing 
the temperature of the ISM (G). 



We also discuss two cases (Sa, Sn) with parameters 
whic h m ay be more appropriate for Seyfert galaxies. In Sec- 
tion 3.2 we consider in general how our results scale with 



the input parameters of the system. 

The range of models presented are chosen to illustrate 
the classification which we have derived f or t he flow mor- 
phologies, and which we discuss in Section 



3.4 



In the figures of this section, we present plots of the 
fiow structures of our simulations. The plots are shown at 
the end of the computational run, or sooner if the fiow 
reaches a steady equilibrium within this time. In all cases, 
the structures shown are broadly characteristic of those 
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Table 5. Parameters of models computed. The labels in column 1 are used to refer to these models elsewhere in the text. See Table |l| 
for definitions of symbols. In all the models computed, half the luminosity is in the anisotropic accretion disc cornponent, /disc = 0.5, 
and the stellar cluster core has uniform density, s = 0. The computational grids for each model are listed in Table W 



Label 


BH mass 




— Cluster 




Mass ratio 


ISM temp. 




Accretion Fig. Comments 






mass, M^i g 


radius, rc.pc 


halo, h 


Mh/M,i 




(/Edd ) 














10 


1 


oo 


- 


1 


n/a 





u 


Pure spherical cluster 


A 


2 


10 


2/3 


5 


0.2 


1 


1 





5 


Fiducial model 


B 


1 


10 


2/3 


5 


0.1 


1 


1 





1 




C 


0.3 


10 


2/3 


5 


0.03 


1 


1 





7 


— Fiducial cluster; 


Ao.005 


2 


10 


2/3 


5 


— 


1 


1 


0.005 


1^ 


— Changes in hole 


Ao.25 


2 


10 


2/3 


5 


0.2 


1 


1 


0.25 


14 


— and accretion 


A' 

^0.25 


2 


10 


2/3 


5 




1 


1 


0.25 1.5 


15 




D 


1 


1 


2/3 


5 


1 


1 


1 





S 




E 


1 


1 


2/30 


5 


1 


1 


1 





9 


Relative 


F 


4 


10 


20/3 


5 


0.4 


1 


1 





IL 


importance of 


G 


1 


4 


3 


oo 


0.25 


10 


1 





u 


the cluster 


H 


10 


10 


2/3 


5 


1 


1 


0.7 





12 




Sa 


1 


0.01 


2/30 


5 


100 


1 


0.01 


5 X 10-* 


17 




Sn 


1 


0.01 


2/30 


5 


100 


1 


0.01 





18 





Table 6. Computational results, for the models defined in Table g. Values are long-term averages; mass transfer rates are in Mq yr~^, 
mass content in Mq. Where values they are marked the simulations which had not fully relaxed (the values given are either rough 
averages for intrinsically variable simulations, or final ones when the flow could not be integrated to convergence). 



Label Afinj Afacc /^acc A/„ind log (iwind/ ergs 1) Mgrid Comments 






9.5 














Mt 


No hole 


A 


8.6 








8.6 


43.3 


1.1 X 10* 


Steady outflow forms quickly 


B 


8.6 








8.6 


39.0: 


105 : 


Many explosions 


C 


8.6 














Mt 


Trapped circulation 


Ao.005 


8.6 


3.2 


0.37 


5.4 


43.1 


9.4 X 103 


Inner disc accretion 


Ao.25 


8.6 


8.6 


1 








8.6 X 103 


Large accretion 


A' 

^0.25 


8.6 


3.7 


0.43 


4.8 


43.1 


8.9 X 103 


Outer disc accretion 


D 


0.86 








0.86 


42.1 


3.0 X 103 


Small cluster 


E 


0.86 








2 X 10-3 


38.7 


Mt 


Small cluster, const iik 


F 


8.6 








8.6 


43.2 


5.2x10^ 




G 


3.8 








3.7: 


43.4: 


1.2 X 10*: 


Hot flow — jet 


H 


8.6 








8.6 





2.2x10* 


Narrow cone — axis force as B 


Sa 


8.6x10-3 


8.6x10-3 


1 








0.30 


'Seyfert' parameters 


Sn 


8.6x10-3 














Mt 


'Seyfert' parameters, no accn 



found throughout the evolution of the flow. The physical 
time which this represents varies, as the length of the run 
was specified in most cases by the total number of compu- 
tational steps taken. 



black hole in the absence of a cluster. As the mass of the 
hole decreases, its region of influence on the global flow con- 
tracts, and the region in which the forces due to the stellar 
cluster dominate consequently expands. 



6.1 Varying the black hole mass 

In the first three models (A, B, C; Figs. ^, ^ and ^, respec- 
tively), we vary the mass of the central black hole from 20 
per cent down to 3 per cent of the mass of the cluster, keep- 
ing the cluster properties and Eddington ratios constant (see 
Table Later, we present the structure of flows in systems 
where the black hole is relatively more massive; however, we 
focus first on this range of parameters since this is where 
mass budget requirements are best met. Also, the flows are 
distinct from the spherical accretion in a massive cluster il- 
lustrated in the previous section and from flows around a 



Model A (Fig. |5|) has the highest mass black hole of the 
first three models, 2 x 10* Mq. The flow relaxes, on a sound 
crossing timescale, to a steady bipolar outflow (cf. Fig. p^, 
which illustrates the relaxation process revealed by the mass 
loss rate from the grid). The inflow in equatorial regions in- 
teracts with the outflow, and is funnelled by a shock into 
a dense region close to the plane. This shock is roughly co- 
incident with the outermost section of the n = 10*' cm-3 
contour (which hits the r axis at 0.2 pc) in Fig. and, in 
Fig. |d, is most easily seen by the sudden change in the di- 
rection of velocity vectors. A narrow, tulip shaped region of 
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a) 



b) 



Danaity/cm^ 



TfllMity/O.Olt: 




r/pc 



r/pc 



Figure 5. Model A at t = 2.6 X 10'* yr. The panels show a) density (from 323 to 1.5 X 10**^001"^) b) total velocity (greyscale and 
contours, up to 1.3), c), d) Magnifications of the central region. The vectors are in the flow direction, and have lengths proportional to 
the velocity, sampled at one vector in every 15 X 15 cells in b), every 5x5 cells in d). The flow is sonic at the u = 0.1 contour. 



circulating, subsonic flow separates the inflow and outflow 
regions. 

The sonic surface separates flow above and below the 
sound speed. In general the streamlines are almost in the 
plane of the sonic surface, although both inflow and out- 
flow from the region are seen close to its ends. Even at the 
ends, the flow is reasonably well-resolved, giving us confi- 
dence that the transition from shock to sub-sonic to super- 
sonic transition through the parallel sonic flow observed here 
is a real feature. Note, however, that this sonic transition 
coincides with the outer radius of the region in which the 
forces are smoothed, and will move inwards if the physical 
'smoothing region' is smaller than that assumed here. 

At the centre of the grid, much of the outflow is con- 



centrated in the grid cells adjacent to the axis (the individ- 
ual cells are 5 x 10"'^ pc square). The flow is confined by 
the inward radial velocities found in these cells (not visible 
in the selected cells shown in Fig. ^ or d), a well-known 
difficulty for axisymmetric hydrodynamical codes. However, 
since the flow here is more-or-less ballistic, and the gravita- 
tional and radiation forces acting on the gas are chosen to 
vary smoothly with angle, this should not have a strong in- 
fluence on the results. This is confirmed by the smoothness 
of the velocity field shown by the contours in Fig. 

Model B (Fig. |) has a central black hole of 10** M©, half 
the mass assumed in Model A. The gravity of the cluster 
and the friction of the mass loading are now sufficient to 
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a) 



b) 



Danaity/cm^ 



TfllMity/O.Olt: 




Figure 6. Model B at t = 6 X 10* yr. The panels are as in Figure pi except that the density range is from 1800 to 2 X 10 cm ^ and 
velocities reach 0.01c. 



stagnate the flow within the cluster. Unsteady structures 
develop, and a sequence of explosions drive out from the re- 
gion. Strong shocks are distributed throughout the region 
where the inflow from the outer cluster interacts with the 
central core; in the plane of the disc (over the surface of 
which a thin plate of very dense, infalling gas forms at 
radii between 0.5 and 1.5 pc), and (as bowshocks) around 
clumps which form in the flow [on axis at z = 1 pc, and at 
{r,z) = (0.12,0.32)]. Where shocks cross each other, regions 
of particularly enhanced density (and pressure) are formed, 
although this pressurization will be transient and may not 
be important for generating cool gas. 

In some cases, these clumps develop an independent 
identity, and fall from the strong outflows close to the axis, 
where they are confined by ram-pressure, towards the plane 



of the accretion disc. However, as they leave the central out- 
flow region the ram pressure of the flow decreases rapidly, 
so the clumps explode and drive shocks into the accreting 
flow. 

We discuss the mass budget of this model further in 
Section 



7.5 



Model C (Fig. ^) has a yet lower mass hole, 3 X 10^ M0. 
The accreted gas collects in a mildly unsteady, convective 
structure at the centre of the grid. Essentially no mass is 
lost from the central parsec-scale region over the 10^ yr for 
which the simulation was run, although the density and size 
of the central structure gradually increases with time. At 
the end of the simulation, the mass contained in the grid 
was 9 X 10^ M0 , and the mean mass loss rate was roughly 
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a) 

Danaity/cm' 




D.S 1 1.5 

r/pe 



c) 

Density/cm 




0.1 0.3 0.3 0.4 0.5 

r/pc 



Figure 7. Model C at t = lO'' yr. The panels are as in Fig. y, exc 
reach 4.6 X 10~^c. 

10~^ M0 yr~^. This total is entirely dominated by gas with 
less than escape velocity flowing out from the edges of the 
accreting flow. Later, once outflows begin to be driven by 
the central convective structure this mass loss rate will begin 
to increase. Even if this were to occur as early as the time 
illustrated in Fig. ^ and if thereafter the mass loss rate were 
to increase in proportion to the flow density (hence roughly 
linearly in time), it would still require 5000 times the length 
of the current simulation before a balanced mass budget 
could be reached. However, it seems likely that the flow will 
remain unsteady in this eventual state, as in Model B, with 
episodic loss of a tiny fraction of the ISM mass. 

The structure of the flow is broadly similar to that of 
Model B. Features to note are the dense plate of gas on axis 
at z — 0.3 pc, which forms where the inflowing halo gas and 



b) 

TfllMity/O.Olt: 




0.1 0.2 0.3 0.4 0.5 

r/pc 



that the density range is from 2330 to 1 X 10^^ cm and velocities 



outflowing wind meet. Oscillations of this plate drive weak 
shocks into the halo, and modulate the rate at which gas 
falls back towards the plane of the accretion disc. The three- 
limbed structure of shocks and sonic surfaces seen here is 
worth noting. It also appears in many other cases (but com- 
pare the single 'tulip' in Model A). Gas injected close to the 
stagnation point does not accelerate monotonically: rather, 
it spirals around in increasing spirals, passing through a se- 
quence of shocks until flnally it can escape. 

6.2 Varying the stellar cluster properties 

We now discuss a selection of models which investigate vari- 
ations in the parameters of the stellar cluster. In Models D 
and E, the cluster mass is reduced. In Model D, we choose 
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the radius of the cluster to be the same as Model B. In 
Model E, we choose the dynamical velocity of the stars to 
be roughly constant. In Models D and F, we investigate the 
behaviour of the flow as the Keplerian velocity decreases 
relative to the sound speed, where F has a rather higher 
velocity of ejection from the central region than D. 

Models D and F can be paired with models A and B, 
respectively, except that they have relatively small stellar 
velocity dispersions (scaled to the sound speed of the flow, 
cf. Fig. I). 

Model G has an cluster with a Keplerian velocity rather 
smaller than the sound speed in the nISM. As a result, the 
flow inside the cluster core is quasi-hydrostatic, rather than 
the supersonic accretion flow found in most of the other sim- 
ulations presented here. Much of the mass loss from the clus- 
ter occurs in a spherical wind. A strong outflow is, however, 
still driven from the centre of the grid, and a recoUimating 
jet is formed as this interacts with the core. 

Model D (Fig. ^ has a lower cluster mass than Model 
B, but the same cluster core radius. Deceleration by gravity 
and mass loading have both decreased, and the flow can now 
drive freely from the centre along the axis of the grid. 

Compared to Model A, the nISM has been made dy- 
namically 'hotter' - i.e. the ratio vk/cs has been decreased. 
The opening angle of the flow which escapes the cluster core 
has been substantially increased. As expected, the sonic re- 
gion substantially expanded; it provides signiflcant support 
against gravity for the gas falling in the accretion disc plane, 
allowing the central outflow to widen. However, it should 
also be noted that with A-L = ALi, the ra,diat,ion forces re- 



tain a Significant influence out to the edge of the cluster. 
This encourages the free outflow of gas along the axis. 

Model E (Fig. ^ also has a lower cluster mass than A, 
and has a far smaller cluster core. Weak shocks are sent into 
the circumnuclear ISM. At the end of the simulation, the 
meiss contained in the grid was 7.6 x 10'^ Mq , and the mean 
mass loss rate was roughly 2 x 10""^ M0 yr~^. This mass loss 
was dominated by gas with sufficient energy to escape the 
nucleus. In this regime, the mass loss rate roughly scales with 
the mass within the grid: a balanced mass budget would be 
expected after 4 x 10® yr, 500 times later than the picture 
presented in the figure. 

Here again we see a three-limbed sonic surface spread- 
ing out through the cluster (note the u/Q.Olc = 0.1 contour 
in Fig. ^ji, which is in this case well within the extended 
smoothing region). At the very centre of the grid, the disc 
infiow shocks into a compact, dense and largely subsonic 
structurp The end nf one nf the three limbs seems tn art as 



rameters to avoid this constraint. However, we have chosen 
here to keep the parameters of the BH and accretion disc the 
same as in Model B, and as a result sacrificed this freedom. 

Model F (Fig. ^ has a cluster with a far larger core radius 
than Model A, but with the same cluster mass. It extends 
our coverage of parameter space to higher ejection velocities 
than Model D, but with a lower cluster dispersion than A. 
Hence, the cluster takes up a far larger fraction of the Parker 
radius of 25 pc. Comparing Figs. |l^ and Fig. ^ the subsonic 
region has inflated from a narrow tulip to take up a broad 
region of the flow. The decreased flow velocities in this region 
have removed the weak shocks close to the accretion disc. 
Funnelling of the flow by the broadened subsonic region has 
acted to offset the decrease in mass flux in the inflow region, 
keeping the density here roughly constant. The velocity and 
opening angle of the outflow are fairly similar between the 
models, so the densities (at a given fraction of the cluster 
core radius) are about 100 times smaller in Model F than in 
Model A - in Model F, the axial outflow is less dense than 
the surrounding cluster ISM, in contrast to most of the other 
models presented here. 

Model G (Fig. 0) has an extremely low cluster core den- 
sity and a very hot ISM. Looking at Fig. |^, we see that Model 
G is an extreme case in our modelling. The Parker radius 
is 1 pc, within the cluster core. If there were no black hole, 
the flow would relax on a dynamical timescale to a subsonic 
core with a sonic transiti on on its surface, hardly a ltered by 
the effects of g ravity (cf. []hevalier & Clegg 1985; Williams 



Dyson 1994). Even when a black hole is present, the sonic 



transition coincides with the edge of the cluster over much of 
its surface, and a substantial fraction of the mass loss flows 
out in this global wind. In Model G, the central outflow is 
highly supersonic. Interestingly, this weakly coUimated cen- 
tral flow is soon recoUimated by the hydrostatic pressure of 
the core: it recoUimates three times within the core, to gen- 
erate a supersonic outflow region with a high aspect ratio - 
a genuine jet. 

Structures similar to these have previously been studied 
in work following on from that of Blandford & Rees (1974; 
see also Wiita 1978a). These authors modelled the outflow 



driven by a central injection of energy, accelerated to super- 
sonic velocities by de Laval nozzles conflned by the pressure 
of the surrounding nISM. The original models treated the 
funnel as very long and thin, but these long funnels were 
found to be unstable to Kelvin-Helmholtz instabilities, re- 
sulting in the lo ss of gas in a series of bubbles rather than a 



a de Lalval nozzle, at (0.004, 0.008), where the flow becomes 



continuous jet (|Gull & Northover 1973; Wiita 1978b 



Nor- 



man et al. 19811). For rather higher injection l uminosities. 



giipprgrinif Ting flnw in anH arniinH tliig ni-.77lp ig irpry ^rarr- 

able. Ttiis results, for the mo.st part, from the variation of 



continuous outflows do occur (Smith et al. 1981 



Smith et al 



the inflow on the disc plane: flow instabilities are expected 
to be suppressed in squat nozzles such as this one dSmithl 
et al. 1S|)83|) . The small mass of the cluster in this simulation 
means that the flow continues to be accelerated beyond the 
nozzle, until it reaches the termination shock which oscil- 
lates about 0.06 pc. 

Note that for the parameters assumed in this case, the 
stellar cluster will be highly collisional. As we have dis- 
cussed, the solution can simply be scaled to different pa- 



198J ). However, in the latter cases the neck of the funnel was 
close to the termination shock of the central wind. In con- 
sequence, the 'jets' so formed were only slightly supersonic, 
and poorly coUimated. Here we find that a dense hydro- 
static envelope can indeed lead to the formation of a well- 
coUimated jet with a substantial Mach number, if a highly 
supersonic but ill-coUimated fiow is driven from the centre. 
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Figure 8. Model D at 6 X 10* yr (steady equilibrium solution). The panels show a) density (from 1513 to 4 X 10* cm~^) b) total velocity 
(greyscale and contours, up to 0.95), c), d) Magnifications of the central region. The vectors are in the flow direction, and have lengths 
proportional to the velocity, sampled at one vector in every 15 X 15 cells. The flow is sonic at the u = 0.1 contour. 



6.3 Varying the distribution of radiation 

Model H (Fig. ^) has parameters identical to those of 
Model B, except for the Eddington ratio. The net force is 
the same along the axis of the disc as in Model B, but is 
outwards in a cone with far smaller opening angle, and the 
inward force in the plane of the accretion disc is far greater. 
Model H quickly relaxes to a steady outflow, albeit with no 
mass leaving the cluster with escape velocity, whereas Model 
B underwent a series of explosions. However, the flow struc- 
ture of Model B close to its long-lived mass minimum of 
3 X lO"* Mq at 2 X 10^ yr (see Fig. is almost indistin- 

guishable from that of the equilibrium state of Model H. 
Altogether, comparison of models B and H indicates that 



the opening angle of the central outflow has little influence 
on the overall flow structure, if the cluster gravity and mass 
input arc significant factors. 

6.4 Accretion over part of the disc 

We next investigate the effect of accretion on the global 
flows. We assume that the accretion from the flow onto 
the thin disc occurs either in its innermost region (Model 
Ao.oos), over a large fraction of its inner surface (Model 
A0.25) or over a large fraction of its outer disc (Model Aq 25). 
These models cover a range of ways in which mass might be 
lost from the global flow to the accretion disc. As accretion 
tends to weaken the net outflow from the central regions. 
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Figure 9. Model E at 9 X 10^ yr (note that this simulation is on a small scale, so the dynamical timescale is short). The smoothing scale 
here is 0.05 pc, the same as in models C, A and B, and takes up a large fraction of the grid. The panels show a) density (from 2100 to 
10^^ cm~^) b) total velocity (greyscale and contours, up to 0.64), c), d) Magnifications of the central region. The vectors are in the flow 
direction, and have lengths proportional to the velocity, sampled at one vector in every 15 X 15 cells. The flow is sonic at the u = 0.1 
contour. 



we use parameters which are perturbations of those used in 
Model A above, which has the strongest initial outflow. 

Models Ao.005 and Aq 25, Figs. ^ and evolve in a 
manner nearly identical to model A, even though around 40 
per cent of the input mass flux is diverted into the accretion 
disc rather than the wind. 

There is a striking similarity between the flows iUus- 
trated in Figs. and [l^. This extends even to the evolution 
of detailed time-dependent features, such as mass clumps on 
axis (Fig. [l6| ) . In one, mass is removed from the simulation 
at the very centre of the grid, in the other across a wide area 
of the outer disc. Some signature of this can be see in the re- 



gions closest to the r-axis in the flgures, beyond r = 0.25 pc. 
Overall, the comparison illustrates convincingly the inde- 
pendence of the form of the solutions from the detailed dis- 
tribution of mass input, or in this case mass removal, so long 
as the total mass budget is similar. 

With its substantially larger area of accretion, the out- 
flow in model A0.25 has been weakened so much that the 
flow structure is more reminiscent of Model C, except with 
dramatically smaller densities in the flow (10^ cm~^ is char- 
acteristic, compared to 10^ cm~^ in Model C), and with a 
far narrower collar of shocked gas at the base of the central 
outflow. 
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Figure 10. Model F at 4.2 X 10^ yr. The panels show a) density (from 350 to 7.5 X lO'^ cm~^) b) total velocity (greyscale and contours, up 
to 1.46), c), d) Magnifications of the central region. The vectors are in the flow direction, and have lengths proportional to the velocity, 
sampled at one vector in every 15 X 15 cells. The flow is sonic at the u = 0.1 contour. 



6.5 Low Eddington ratio 

The models described so far in this section were for bright 
QSOs. We now consider the effects of changes in parameters 
to values characteristics of Seyfert nuclei. 

The similarity in line widths and profiles between AGN 
of widely differing luminosities can be interpreted as im- 
plying that thp Vrripmatrna arp gimrlar Tf thig ig in 

Seyfert Inuclei the characteristic length and timescales must 
be A ~ r ~ 0.01, relative to the QSO case. If the Edding- 
ton ratio, etc., remain constant then such a self-consistent 
Seyfert model would have masses smaller than those of QSOs 
by A, and densities increased by 1/A. The ionization param- 
eter remains constant so long as Q is determined to within 



a small factor by stellar physics. In this case, the ratio be- 
tween cooling and flow times is constant, and the models 
scale directly from those presented above (see also the next 
subsection) . 

However, many authors find black holes masses in 
Seyferts ~ 10* Mq (as might be expected in an old QSO 
brought ba ck to life), requiring that a ccretion is sub- 



Eddington ( Padovani, Burg fc Edelson 1990 



Lasota et al 



1996 ). Structures similar to those shown in the preceding 
subsections would persist if the efficien cy o f outward radia- 
tion driving were increased (cf. Section 7.4). 

However, if there is no such enhanced driving the topol- 
ogy of the flows changes radica lly, as we have discussed from 
scaling arguments in Section 3.3.2 In this subsection, we 
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Figure 11. Model G at 9 X 10 yr. As this simulation was very time-consuming, it was halted at this stage, before the solution had 
reached equilibrium. However, since the wind mass loss has reached 97 per cent of the injection rate, and little movement was seen in 
the solution, these plots should be fairly representative of the eventual equilibrium. The panels show a) density (from 211 to lO'^ cm""^) 
b) total velocity (greyscale and contours, up to 2.7, note that for the gas temperature of 10* K used here, the sound speed is 0.32 so this 
corresponds to a peak Mach number of 8.7), c), d) Magnifications of the central region. The vectors are in the flow direction, and have 
lengths proportional to the velocity, sampled at one vector in every 32 X 32 cells in b), every 7 X 7 in d). A jet close to the axis passes 
through recoUimation shocks which reach the axis at 1, 2 and 3pc before it escapes the cluster core. 



present numerical results for the case in which there is no 
increased driving. We consider a black hole mass of 10* M©, 
but a cluster mass of only 10® Mq , sufficient to power ac- 
cretion at luminosity lO^^ergs"^, i.e. about O.OlI/Edd- (We 
shall see that the form of the flow is similar to that for a 
low mass-loss cluster of 10* M0 , as might be expected in 
a 'fossil QSO' nucleus.) In order that stellar collisions do 
not dominate the mass input, the core radius of the clus- 
ter must be larger than 0.07 pc (note that in this case the 
stellar velocities will be dominated by the mass of the black 



hole) . This core radius is an order of magnitude greater than 
variability-based estimates of the size of the BELR (e.g. Pe- 
terson 1994), which suggests that the BELR lies entirely 
within the cluster core. 

We calculated two models, with and without central 
accretion (Figs. ^ and |l^, respectively). In both, we found 
near-spherical inflows of the nISM throughout the nucleus, 
except for a small hydrostatic core in model Sn. 

These low Eddington ratio models do not show the 
bipolar structures characteristic of the QSO models, and 
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Figure 12. Model H at 2 X 10* yr. The panels show a) density (from 400 to 2 X 10^^ cm~^) b) total velocity (greyscale and contours, up 
to 3.3), c), d) Magnifications of the central region. The vectors are in the flow direction, and have lengths proportional to the velocity, 
sampled at one vector in every 31 X 31 cells. The flow is sonic at the u = 0.1 contour. 



as observed in radio maps of Seyfert galaxies. This is be- 
cause the accretion disc luminosity is only a fraction of the 
(insignificant) radiation force. Spherical inflow will charac- 
terise Seyferts with /Edd <C 1 if the opacity is dominated by 
Thomson scattering. Observations imply that other sources 
of opacity will produce bipolar flows in Seyferts. The signa- 
ture of this radial inflow should be detectable in the vari- 
ability of line profiles, with the red wings of the lines leading 
the blue. For the smaller nuclei observed in Seyfert galaxies, 
indivi dual supernovae will often evac uate the entire nuclear 
ISM (Perry, Williams & Dyson 1999) returning the fiow to 



initial conditions similar to those we have assumed for our 
simulations. 



7 DISCUSSION: THE PROPERTIES OF 
FLOWS IN ACTIVE NUCLEI 

We have modelled the central region of an AGN, assuming 
only that a black hole and accretion disc with an aspheri- 
cal radiation field is embedded in a young starburst stellar 
cluster. These two simple components give rise, symbioti- 
cally, to a nuclear interstellar medium with a remarkable 
and fascinating variety of structures. Despite the absence, 
in our models, of any effects due to radio-related phenom- 
ena, the nISM shows a range of behaviour normally asso- 
ciated with relativistic effects near the black hole: episodic, 
explosive, percolating structures, and/or well coUimated, hy- 
personic hydrodynamic jets - given only that the radiation 
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Figure 13. Model Aq.oos at 1.8 X lO** yr (steady equilibrium solution). The panels show a) density (from 191 to 1.1 X lO^^cm"^) b) 
total velocity (greyscale and contours, up to 1.3), c), d) Magnifications of the central region. The vectors are in the flow direction, and 
have lengths proportional to the velocity, sampled at one vector in every 15 X 15 cells. The flow is sonic at the u = 0.1 contour. 



field drives an outward flow along the axis of the disc. When 
the Keplerian velocities of the stars, vk, are ^ Cs, the black 
hole drives a meridional flow within the previous hydrostatic 
core. The mass-loaded outflow from the core is either con- 
fined by an (apparently very unstable) internal termination 
shock or can drive gas beyond the cluster core, in episodic 
explosions or as a continuous wind. For smaller vk/cs, con- 
fined jets can form on-axis. In all cases where a wind drives 
out through the nucleus, it has a half-opening angle which 
is rather smaller than that of the (conical) region close to 
the centre within which the net force is outwards (typically 
30 degrees compared to 60 degrees). 

The models we present here are distinguished from pre- 
vious studies by the effects of the stellar cluster and by the 



nature of the flow at small radii. The poorly coUimated cen- 
tral wind from the very centre of the grid drives outwards 
supersonically into the cluster. When this outflow has suffi- 
ciently great momentum, it blows straight through the clus- 
ter and onwards into the galaxy ISM: such fiows soon settle 
to steady equilibria. For lower momentum flows, the out- 
flow shocks against the inflowing nISM of the cluster, and 
falls back to the core in a circulation. In some cases, a series 
of explosions are found; rather than buoyancy-driven bub- 
bles, these are pushed outwards by the recycling of matter 
through the black hole core. The timescale for these explo- 
sions is the dynamical timescale of the nuclei, which will 
vary between roughly 300 yr in QSOs and 3yr in Seyfert 
galaxies. These times, too long to explain continuum vari- 



© 0000 RAS, MNRAS 000, 



30 R.J.R. WilUams, A.C. Baker & J.J. Perry 



a) 



b) 



Danaity/cm^ 



TfllMity/O.Olt: 




r/pc 



r/pc 



Figure 14. Model A0.25 at 5.1 X 10* yr. The flow here is oscillating, cf. Fig. The panels show a) density (from 171 to 5.1 X 10* cm~^) 
b) total velocity (greyscale and contours, up to 1.2), c), d) Magnifications of the central region. The vectors are in the flow direction, 
and have lengths proportional to the velocity, sampled at one vector in every 15 X 15 cells. The flow is sonic at the u = 0.1 contour. 



ability, may relate t o variations in the underlying structure 
of line-emitting gas (Perry, van Groningen fc Wanders 1994 



Wanders & Peterson 199(: ) 



Gas which rains onto the disc in many of our models 
may also suppress the wind from some part of its surface. 
If the gas in this region can cool, it will irradiate the outer 
disc and may have important ramifications for angular mo- 
mentum transport in the disc. 

The structures we predict for the flows are very sensi- 
tive to the central luminosity, if it is close to the (on-axis) 
Eddington limit. The boundary between stagnated inflow 
and free, steady outflow is fine, depending on whether tran- 
sient explosions can drive significant mass from the nucleus. 
Indeed, there seem to be bimodal solutions in some cases. 



where the fiow can take either a low-mass, free-outfiow or 
high-mass, recirculating form. Quite small changes in lu- 
minosity may lead to ejection of large masses of nISM gas 
within the dynamical timescale of the nuclear fiow. In a fully 
self-consistent model, variations may be even more catas- 
trophic than those presented here (although these variations 
may be damped by lags due to viscous transport within the 
accretion disc). 

In many cases the time-steady system does not lead to 
structures in the nISM which are time steady, unlike the 
models of the wind from the surface of a molecular torus 



by Balsara & Krolik ( 1993 ). Our results indicate that the 
fiows may remain chaotic for all time, with shocks forming 
an important part of the structure. In the models of Balsara 
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Figure 15. Model Aq 25 at 4.6 X 10* yr. The flow here is oscillating, cf. Fig. The panels show a) density (from 182 to 9.6 X 10^ cm"'') 
b) total velocity (greyscale and contours, up to 1.3), c), d) Magnifications of the central region. The vectors are in the fiow direction, 
and have lengths proportional to the velocity, sampled at one vector in every 15 X 15 cells. The flow is sonic at the u = 0.1 contour. 



& Krolik, outflows are driven from scales ^ 1 pc - the flow 
within the BELR is purely accretion - and have significantly 
lower velocities than those we find in our models. 

This is also true of t he Co mpton-heated disc winds mod- 
elled by Woods et al. (1996). In their very detailed treat- 
ment, they found low- level variability in the flows - it is our 
eventual aim to model the thermal structure of the flow in 
similar detail, but cooling of gas behind global shocks will 
make this problem even more computationally taxing than 
the cases calculated by Woods et al. 

We now consider the physical properties of the nISM 
flows found above. We first discuss the free parameters which 
are most important in determining the flow structures, and 
the physical limits of the models. Next, we derived the dy- 



namic pressures and ionization state of the flowing gas, espe- 
cially in Section 7.3, where we also consider the interaction 



of the nISM with the shocks driven by stars and supernovae. 
We calculate the overall mass budget of the nucleus in Sec- 
tion 7.5 and the mechanical energy carried in the winds and 
their effects in Section 7.6. Throughout this section, we high- 
light the observational consequences of the model. 



T.l Classification Scheme 

Comparison of the fiow structures shown in the previous 
section with the sketched structures in Fig. ^ illustrates our 
suggested classification in terms of the first two parame- 
ters. In terms of the physical parameters of the nuclei, the 
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Figure 16. Comparison of (a, b) model Aq.oos and (c,d) Aq 25, both at 6.5 X 10^ yr. Notice the similarity of detail in these simulations 
(away from the disc plane), in particular the dense clump of gas at 2: ~ 0.65 pc on axis. 



Keplerian velocity is = GMc/rc and the sound speed is 
determined by the shape of the AGN continuum. We can es- 
timate the ejection velocity, v^j , by assuming that the central 
outflow is dominated by gas accelerated in the z direction 
within the central smoothing region, so 



^cj _ r n t 1 Gh 'h \ 



1/2 



(21) 



The comparison between the calculated and the measured 
Wej is shown in Fig. ^. The comparison is good except for the 
cases (described in the figure caption) where the approxima- 
tions made in applying eq. ( pl| ) break down. 

Where the outward driving from the black hole is large, 
the flow relaxes to a steady bipolar form: gas falls inwards 
in the plane of the disc and is then driven away when 



reaches the super-Eddington region on-axis, close to the cen- 
tral black hole. When the retarding effects of cluster gravity 
and mass loading become great enough, the outfiow shocks 
within the stellar cluster. If the sound speed in the hot gas is 
sufficiently great, a steady recoUimating jet structure can be 
produced. More characteristic, however, is a structure where 
episodic explosions drive gas from the nucleus in bursts. As 
the retarding forces increase further, the outffow region be- 
comes more compact and can less often drive mass out from 
the nucleus, eventually becoming a turbulent core in the 
flow onto which gas injected in the outer regions of the stel- 
lar cluster accretes. 

For our approximations to hold, we require that the 
forces due to the black hole dominate close to it, but that 
the gravity of the cluster dominates at large radii (this as- 
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Figure 17. Model Sa at 2.5 X 10 yr. Notice the smaller scale of the grid in this simulation: the flow reaches equilibrium as early as 200 yr 
as a result of the very short dynamical timescale in this model. The panels show a) density (from 15 to 2 X 10* cm~^) b) total velocity 
(greyscale and contours, up to 1.2), c), d) Magnifications of the central region. The vectors are in the flow direction, and have lengths 
proportional to the velocity, sampled at one vector in every 15 X 15 cells in b), every 5 X 5 in d). Gas flows freely in from the cluster 



and is accreted at the centre of the grid. The anisotropic radiation, 
perturbations on this flow. 



with a mean Eddington ratio of only 0.01, produces no noticeable 



sumption is strained in models D and E). Various other pa- 
rameters come into play; although they are less important, 
they limit the strict application of our classification scheme. 
Changes in the opening angle of the outflow (Model H) and 
the distribution of stars within the cluster will have only a 
small effect. Accretion of gas from the flow on a dynami- 
cal timescale can also serve to weaken the central outflow 
(models Ai). In all these cases, the morphology of the flows 
are still broadly within our classification, and analytic es- 
timates of fiow velocities and densities are still readily cal- 
culable. The net effect of accretion on the fiow structures is 



to weaken the winds and dampen the variability, although 
this depends on where the accreted mass is removed from 
the global fiow. Changing this region does not fundamen- 
tally alter the structure of the fiows until the fraction of the 
mass accreted exceeds ~ 0.5. The overall fraction is more 
important in determining the fiow structure than the loca- 
tion of the sink. The structures are basically similar to those 
found in non-accreting cases, albeit with different black hole 
parameters (cases in which there is a central wind blowing 
through a mass-loading cluster will again be similar, so long 
as this wind does not have a high angular momentum). 



© 0000 RAS, MNRAS 000, 



34 R.J.R. Williams, A.C. Baker & J.J. Perry 



a) 



b) 



Danaity/cm^ 



TfllMity/O.Olt: 




O.OZ 0.03 
r/pc 

Figure 18. Model Sn at 4.2 X 10^ yr. Notice the smaller scale of the grid in this simulation: the flow reaches equilibrium as early as 200 yr 
as a result of the very short dynamical timescale in this model. The panels show a) density (from 15 to 1 X lO^Ocm"^) b) total velocity 
(greyscale and contours, up to 1.5), c), d) Magnifications of the central region. The vectors are in the flow direction, and have lengths 
proportional to the velocity, sampled at one vector in every 15 X 15 cells in b), every 5 X 5 in d). Gas flows freely in from the cluster 
and is accreted at the centre of the grid. The anisotropic radiation, with a mean Eddington ratio of only 0.01, produces no noticeable 
perturbations on this flow. 



This classification should be supplemented by the num- 
ber of dynamical times for which mass is retained by the 
nucleus (which will often relax to a constant value). Most 
importantly, the gas may not remain in Compton equilib- 
rium, in which case the scalings we have derived will have 
to be recalibrated in terms of an 'effective' Eddington ratio. 

7.2 Densities and velocities 

The models we have calculated imply a range of local prop- 
erties (densities, pressures, velocities) for the flows. We now 



compare these with simple estimates, based on spherically 
symmetric flows, similar to those often given in the liter- 
ature. Conservation of mass in spherical symmetry gives a 
density for the global flow (on scales of r = rpc pc) of 

n ~ 4 X 10 -5-^ cm (22) 

\Macc77acc,-l / ^Jc («flow/0.01c) 

where u is a characteristic flow velocity within the nucleus 
(by comparison, for the sound speed, v = lQ~^Tj^^c and for 
the Keplerian velocity, « ~ 7 x 10"^(Afci,8/10rpc)^''^c). 
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This estimate of density is in line with the values seen 
in the accretion flows in, e.g., model A. The flow velocity in 
the inflow region is about 3 x 10~^c, similar to the Keple- 
rian velocity and the density in this region corresponds to 
equation to reasonable accuracy. In the axial outflow 
region, the velocity is about 3 times higher (corresponding 
to non-Keplerian velocities that are obviously required by 
the widths of the observed optical emission lines) . Here, the 
density is higher than in the inflow region, with a value more 
than ten times larger than suggested by equation (^) . This 
is because the 60 per cent of the mass input not accreted by 
the black hole flows out in a cone of half opening angle of 
only ~ f5°, (3.5 per cent of the solid angle) rather than the 
large fraction assumed in deriving the approximate result. 

As the outflows are driven from very close to the central 
black hole, their velocities are characteristic of these radii. 
While the highest speeds in the models presented here are 
~ 8700 km s~^, this value can be rescaled by choosing dif- 
ferent flow parameters, and might be significantly increased 
if the effective smoothing length was lower than that as- 
sumed. The physical processes whic h dete rmine the smooth- 
ing length are discussed in Section 2.2.6. It is possible that 
values as small as ~ 3 x 10^^ cm may be appropriate in 
some cases, rather than the values of ~ 1.5 x lO'^^ cm we 



used principally because of computational necessity. Nuclei 
with very tiigii outflow velocities are, however, likely to have 
structures similar to our model A, with little retained gas 
and few global shocks. 

If mass is retained by the nucleus for N dynamical 
times, as for instance in models B and C, the density and 
pressure [equations (^^ and (^^] will be boosted by a factor 
A'^, and the ionization parameter [equation (^H])] decreased 
by a similar factor, within the central convective structure. 
Where this occurs, this will increase the emission from gas 
in the central plume beyond that, retention occurs. 

Therefore, we find that the simple arguments of PD 
and Perry (1993a) are borne out throughout much of the 
nucleus. However, they did not anticipate several important 
effects which we have found boost the gas density in some 
regions. These are the very regions which may dominate 
the line emission from these nuclei (as we have seen in the 
discussion of Fig. M above). 



7.3 Stagnation pressures, ionization parameters 
and the multi-component BELR 

Both the ionization parameter and the stagnation pressures 
within the nISM provide checks on the consistency of our 
models with our assumptions, and are first steps towards 
assessing the observational consequences of the symbiotic 
starburst- black hole model for AGN. In this section, we 
present results for the ionization parameter in the freely 
fiowing nISM, and the ionization parameter and gas pres- 
sure at stagnation in the nISM, for two characteristic shock 
velocities. 

The ionization parameter we use is that as defined by 
(Krolik et al. |l98l| ) 

Fi 



nuksTc' 



(23) 



where _Fi is the local radiation flux between 1 and 1000 Ryd. 
This choice of ionization parameter is that which is most 



appropriate for discussions of the interplay between hydro- 
dynamics and photoionization equilibrium. Collin-Souffrin 
(1993) discusses the relationship between the various ion- 
ization parameters in common use. We plot H based on the 
local intensity of radiation, assuming that the ionizing fiux 
is a constant fraction of the total flux. 

PD and Perry (1993a) emphasized the importance of 
shocks in the flows in AGN. A simple estimate of the pres- 
sure on cool gas behind a stationary shock, i.e. the stagna- 
tion pressure of the flow at Vs, is 



Ps 



4 X 10' 



/e 



/^acc^acc, — 1 



o 



■ cm"'^ K (24) 



and the ionization parameter of shocked gas is 

CUflow 



(25) 



Hot phase gas will begin to cool when its ionization 
parameter falls below that at the knee of the equilibrium 



100 



Li 



T, 



-■ill 



0.1 



(26) 

cooling begins to domina te over 



ifiol 

where bremsstrahlung 
Com ptonization. For a recent composite spectrum, (Ferland 



199e| ) Li/iBoi = 0.45 and Tc.i = 1.3. 

While this limit is consistent with our choice of con- 
stant Compton temperature in the flows we analysed, the 
anisotropic component of the radiation field may well have 
a different spectrum to the isotropic one. For example, for 
the outer regions of the disc plane, it might be more ap- 
propriate to consider only the luminosity under the blue 
bump. In the case of Ferland's continuum, the luminosity in 
the 1-10^ Ryd wave-band then decreases by a factor more 
than 5 with a concomitant decrease in 5 (and, as discussed 
below, the Compton temperature increases by a factor 2). 
These gross changes in spectral shape may have consider- 
ably more far-reaching effects than simply a change in the 
ionization parameter. The effects of variations in spectrum 
and Compton temperature with angle will be discussed in 
detail in future papers. 

If Hf 10, our assumption of Compton equilibrium for 
the hot gas will be strained by the increasing influence of 
bremsstrahlung and atomic cooling. Conservatively, we take 
Hf = 10 as a limit below which cooling is likely to become 
important. 

The pressure in any gas which manages to cool will 
rarely be less than that in the surrounding free flow, and 
so the free-flow ionization parameter Hf is an upper limit to 
that which acts on the cool gas. This free-flow Hf is shown 
in panel (a) of the graphs in this section. 

In panel (b) of the graphs in this section, we show the 
ionization parameter Hg in gas shocked to stagnation in the 
rest frame of the nucleus. This ionization parameter char- 
acterises the 'average' ionization parameter in shocks which 
form around local obstacles to the flow, as well as in steady 
features within the flow. 

In reality, most of the obstacles in the flow will be driven 
by stars and supernovae. The speed of shocks driven by these 
stars will depend on their orbital velocities, mass ejection 
velocities and the amount of gas swept up by the shell, as 
discussed extensively by Perry (1992, 1993a, b, 1999). As an 
indication of the likely importance of these random veloc- 
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_i I ■ ^ ■ 



b) 



stagnation 1 




Figure 19. Model A at 2.6 X 10* yr. The panels show a) ionization parameter, S, in the free flow, b) ionization parameter, Ss, at 
stagnation in the rest frame of the nucleus, c) 'minimum' stagnation ionization parameter, Sj;, at stagnation against a counter-moving 
Keplerian obstacle, d) gas pressure, ps, at stagnation in the rest frame of the nucleus (in cm~'^ K). 



ities, we also plot ionization parameters behind shocks of 
velocity vhow + vk , where vk is the local velocity of a circu- 
lar orbit in the potential of black hole and cluster. This is 
the maximum shock velocity for a rigid obstacle. Note that 
in the inner regions of the flow, where the gravity is dom- 
inated by the black hole, the Keplerian velocities depend 
on the absolute mass of the central black hole (and thus 
the degeneracy between assumed values of (/Edd) and A/h is 
lifted). 



In Fig. |19|a we show the ionization parameter, Ef, for 
Model A. As Hf > 10 through most of this region, our as- 
sumption that the flow is isothermal at the Compton tem- 
perature is self-consistent, except for a very small region 
close to the origin. There are no direct observational tracers 



of gas with this high temperature and ionization param- 
eter. Even iron K-shell absorption, one component of the 
so-called 'warm absorbers', which can be significant in gas 
with temperatures as high as 10^ K, will only be important 
in gas with typical densities ^ 10^ cm""^ for the radiation 
flux assumed here (Mathews & Ferland 1987). 



While a wide range of densities and pressures ar e pre- 
dicted by these models (cf. empirical evidence for this 



Bald- 



win ct al. 1995), broad systematic trends are also expected. 



For instance, in Model A, there are two regions where the 
gas shocked to typical stagnation pressures can form line- 
emitting clouds. In the rest of the flow, post-shock gas has 
high ionization parameters, and so does not cool sufficiently 
to radiate emission lines (Fig. |l|b). The two line-emitting 
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regions are near the outflow cone, and just above the disc. 
The stagnation ionization parameter is ;^ 0.1 through most 
of the conical outflow around the axis. The pressures here are 
high lO'^* cm~^ K, Fig. pj|d). Above the accretion disc, 
in the other low ionization parameter region, inflowing gas 
shocks to a rather smaller pressure, ~ lO'^^ cm"'' K. The rela- 
tive contributions of these two regions to the observed broad 
emission lines will be weighted by their relative volumes of 
rotation about the z-axis. Many physical effects other than 
the overall pressure and ionization parameter must be in- 
cluded in determining the volume emissivity of the cool gas 
from hot phase properties (cf. PD, and below). However, the 
obvious distinction between the properties of gas in the two 



regions means that, broadly-defined, two 'components' will 
describe these two, dynamically separate flows. 

The fastest shocks, where Keplerian velocities add to 
the flow velocity, can cool gas throughout the entire cluster 
core as well as in a broad cone outwards along the flow axis, 
if the gas can be maintained at the high pressure at the 
head of the sh ock for at least a cooling time (Fig. pj|c). In 
future papers (Albertsson & Perry 1999), we will predict the 



detailed line profiles and variability signatures expected as 
a result of the combination of all these effects. 

The flow in Model B (Fig. |2c|a) is too dense in most of 
the central 'plume' to be maintained in isothermal equilib- 
rium by Compton heating and cooling. Except for a small 
conical area about the central black hole, all the gas with 
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Figure 21. Model Sn at 4.2 X 10^ yr. The panels of this plot show a) ionization parameter, S, in the free flow, b) ionization parameter, 
Hs, at stagnation in the rest frame of the nucleus, c) 'maximum' stagnation ionization parameter, Sjj, d) gas pressure, ps, at stagnation 
in the rest frame of the nucleus (in cm~'^ K). The anisotropic form of the ionization parameter plots results from the variation of nuclear 
luminosity and spectrum with angle to the nucleus. 



r ^ 0.5 pc and z 1 pc has an ionization parameter less 
than unity. This low free-flow ionization parameter puts 
considerable strains on our assumption of Compton equilib- 
rium. Net heating of the gas in smaller mass-loading shocks 
may be able to offset some of the bremsstrahlung cooling. 
In addition, the additional cooling will have the side-effect 
of increasing the opacity of the ga s, a nd hence the radia- 
tive forces on the flow (cf. Section 7.4). As these radiative 
forces increase, the topology will gradually become more like 
Model A. This self-regulation will enhance the frequency 
with which flow structures intermediate between Models A 
and B are found. The particularly high pressures found in 
the accretion disc plane will encourage cooling, perhaps lead- 



ing to the addition of mass from the free flow to the outer 
accretion disc, as in Model A0.25. 

As might be expected, the distribution of stagnation 
ionization parameters, in Model B (Fig. pc|b) is rather 
smoother than the distribution of free-flow ionization pa- 
rameters discussed above (since gas at rest has the stag- 
nation pressure by deflnition, slow moving shocks have a 
less marked effect on the ionization parameter). The stag- 
nation ionization parameter is, however, less than 1 through- 
out the central 1 pc sphere of the nucleus. In the near-radial 
inflow region outside the central plume, gas pressures are 
^ 10^^ cm-^K. 

The stagnation ionization parameters are far lower than 



© 0000 RAS, MNRAS 000, 



Hydrodynamics of the ISM in Symbiotic AGN 39 



this in the central plume itself, with pressures reaching more 
than 10^^ cm~"^ K. Shocks in the plume will cool fast, and 
the cool gas will radiate strongly, since this is the region in 
which the radiation from the central accretion disc is at its 
strongest. 

In Fig. ^ we show the ionization parameters de- 
rived from the low Eddington-ratio model with no accre- 
tion, Sn (cf. Fig. |l^). The free-flow ionization parameter, Ht 
(Fig. pl|a), shows that the nISM is always hot and diffuse. 
The stagnation ionization parameter, Hg (Fig. |2llb), is suffi- 
ciently small to allow cooling in only a very small region in 
the plane of the accretion disc where the radiation field is 
weak. This region is expanded to fill most of the cluster core 
when the infiuence of Keplerian velocities is accounted for 
(Fig. |2l|c) - however, the cooling length behind such shocks 
would be a substantial fraction of the size of these nuclei. 

B ehind shocks in this region the cooling will be very 
rapid (Innes & Perry 1999), and consequently the size of the 
BELR will b e close to that derived from variability studies 
( |Perry 199^ ). 

W e conclude that the BELR in many AGN must be 



highly complex, anisotropic, unsteady and irregular. Physi- 



cal effects not included in the current treatment will modify 
our results, but the general form of the flows will remain 
similar to those presented here. Nevertheless, Figs. ^ and 
give a general impression that the structures should be 
able to be broken down into a few components, with vary- 
ing mean properties. We will undertake the detailed mod- 



elling required to confirm this in future papers (Albertsson 
& PerrI 19991). 



7.4 Column density and opacity 

In some of our models, in which mass is retained for several 
dynamical times, the column depth of the nISM reaches val- 
ues around 10^** cm~^. This will have potentially observable 

pnrpK At t.Tipsp rnlTimn densities the niirlpiiq will 



conseqi 

be gett ing close to a Thomson optical depth of unity. Elec- 
tron scattering will weaken the continuum observed along 
the densest lines of sight to the nucleus, smooth out con- 
tinuum variability and lead to polarisation of the scattered 
component. The emission lines will also develop a thermally 
broadened component, although the optical depth to line- 
emitting clouds distributed in the nISM will be smaller than 
that to the central continuum source. 

X-ray absorption by highly ionized oxygen and iron, at 
energies of 0.3-10 keV, is s een in many low luminosity AGN 
(iNandra fc Pounds 199^ ; [Reynolds fc Fabian 199^ ; [Krolik 



& Kriss 1995). The location of the absorbing gas is uncer- 



tain: some of the absorption may be from the gas which is 
observed in Seyfert 2 galaxies to scatter the obscured cen- 
tral radiation source, while some may also be far closer to 
the nucleus. This 'warm absorption' is currently the best 
diagnostic of gas which is more widely distributed than the 
locally confined broad emission line clouds. 

The observational results are, in general, given in terms 
of the column density and ionization parameter required for 
a uniform intervening slab of solar-abundance gas in thermal 
equilibrium to generate the absorption features seen. Typ- 
ical column densities are cm and temperatures 



lO'^-lO^K. Reynolds & Fabian (1995), for instance, base 



a column density of 10^^ cm~^ and Li/nr^ — 30ergcms~^, 
where Li is the total ionizing luminosity and n is the elec- 
tron density. These values imply that the gas temperature 
is 10^ K - they find the gas resides in a toe-hold of thermal 
stability between the usual cool and hot equilibria. The poor 
spectral resolution of (all but the most recent) X-ray obser- 
vations is a valid reason for using this appro ach to get the 
most out of the available data. Krolik & Kriss (1995) empha- 



sise that the material observed may well be out of thermal 
equilibrium, cooled by adiabatic expansion (although it is 
likely to be in ionization equilibrium). This is an even more 
important constraint where the gas is being evaporated from 
a local evaporation centre rather than expanding as part of 
a global flow. 

Until our models are extended to explicitly include the 
thermal balance of the gas, we cannot model the warm ab- 
sorber in detail. It is, however, clear that a number of poten- 
tial locations for the warm absorber arise naturally in these 
models. The global flow may dense enough that it has begun 
to cool in some regions, or may be cooled a t large distances 



from the nucleus by adiabatic expansion (Krolik & Kris; 
1995). In addition, in the shocks which create the broad- 
line emitting gas must cool through this temperature range; 
subsequently the gas will evaporate from the cooled clouds 
and rejoin the flow. For a plane radiative shock, the col- 
umn between lO'' K and lO'^ K is roughly equal to the mass 
flux multiplied by the cooling time, which is J5 10^^ cm~^ 
for post-shock pressure of 10^* cm~'^ K. The Compton heat- 
ing time is ~ 10^ times longer than this cooling time, and 
may be lengthened further by adiabatic cooling and resid- 
ual atomic cooling. However, the evaporation will only occur 
once the gas has reached a region with substantially lower 
pressure than behind the radiative shock and further diver- 
gence of streamlines will occur while the gas is heating (PD), 
both of which effects may act to decrease the column of gas 
available to generate absorption. The he ating gas may be a 



signi ficant component of the absorption ( Reynolds fc Fabian 



1995), but this needs to be studied in greater detail. 



7.4-1 Cooling and the effective Eddington ratio 

Cooling leads to significantly enhanced opacity and hence 
radiative forces on nISM gas: in consequence, the densities 
will be reduced, so an equilibrium may be found with a 
substantial mass loss rate even where the Eddington ratio is 
less than unity. The overall radiation force can be expressed 
as an angle-averaged effective Eddington ratio, 

.?7accMacc ^ „ Md 

where r^acc is the efficiency with which the central accre- 
tion disc radiates the rest-mass energy delivered to it by the 
global fiow; /iacc is the fraction of the mass input from the 
cluster delivered to this central region, and a is the mean 
effective cross-section for gas close to the black hole. Of all 
the various uncertainties implicit in the terms of the above 
equation, the determinat ion of the effective cross-section is 
perhaps the greatest (cf. 



/ /■Effcctivc\ ^ / -f \ 
UEdd / — \/Edd) 



Turner et al. 1993| ) 



their discussion of the observations of a particular Seyfert on 



In most of the models presented here, we assume that 
the overall Eddington ratio is close to one. Various indepen- 
dent observations suggest that typically 10~^ < fF,dd < 1, 
given specific model assumptions (see Appendix 2.2.4). As 
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/Bdd approaches unity, the structure of the accretion disc 
will become bloated by radiation pressure, although this 
bloating is reduced if the luminosity of the accretion disc 
is dominated by reprocessing. 

We have assumed ct = ctt (and thus /edd = /fidd, 
the classical Eddington ratio) in the models presented here. 
However, if the density of the gas is high enough, some frac- 
tion of the gas may be able to cool, enhancing the radia- 
tive driving. For optically t hin gas at equ ilibrium in the 
BELR (modelled by Cloudy Ferland 1996) the force mul- 
tiplier a-/(TT increases to roughly 2 at lO'' K, 40 at 10^ K 
and to 10^ and beyond when the gas has cooled fully (see 
also lArav, Li & Begelman 1994). It is reasonable to take 



q = 7 X 10^/0 when the gas is cooler 1: han ~ 5 x 10^ K (cf. 



Roser 1979 



Dyson, Falle fc Perry 1981 ) , although this force 



is dominated by resonance line scattering and will weaken 
by a factor of a few if the lines are optically thick. Hence 
even in Seyfert galaxies, where /sdd (determined by elec- 
tron scattering) may be far smaller than unity, a dynamical 
feedback may lead to flows of very similar forms to those 
presented in here. 



7.5 Mass budgets: winds and accretion efficiency 

Here, we discuss the mass budgets found for the models we 
have calculated: how much of the mass lost by the stellar 
cluster - and its associated heavy elements - remains in the 
nISM, how much is delivered to the accretion disc, and how 
much is transferred to the wider galaxy. Numerical values 
for the equilibrium mass processing rates have been given in 
Table ||. 

The mass processing rates of models A and B are shown 
in Fig. ^ (we do not show Model C, since almost no mass 
is lost from this nucleus over our integration period). For 
Model A, the mass loss rate reaches an equilibrium after 
roughly the sound crossing time of the nucleus. Model B 
produces a series of explosive mass loss episodes, continuing 
for times as long as we have simulated. The characteristic 
timescale of these explosions, typically 10"^ yr, is comparable 

o 2 /2 1/2 

to the gravitational timescale, 1.6 x 10 rc'pcMg yr . 

We plot the variation of total mass within the grid of 
Model B in Fig, p^ . This is the integral of the variation 
shown in Fig. p3b. As we have already described, the flow 
is highly variable, alternating between times when the in- 
put mass is almost entirely retained (the steepest upward 
slopes in Fig. p4[ ), and epochs of rapid mass loss. In particu- 
lar, this figure emphasises the long-term chaotic behaviour of 
the flow structure we have simulated. It also shows how close 
Model B is to the border between steady-outflow and chaotic 
flow. During the extended minimum around 2 x 10^ yr, the 
structure and density of the flow is very similar to that 
shown, for Model H, in Fig. ^ Indeed, we would hope this 
were so, since the fundamental parameters for Models B and 
H are identical (cf. Fig. |^ - they differ only in the opening 
angle of the central outflow cone. 

The fraction of stellar mass loss delivered by the flow to 
the accretion disc depends critically both on the structure 
of the flow and on details of the processes by which mass is 
added to the accretion disc. A primary aim of this paper is 
to address the first of these issues. Therefore, we treat the 
second by three extreme assumptions - no mass addition to 



the disc, or addition only in an inner or outer region of its 
surface. 

The fractional rate of accretion, /iacc, was assumed to 
be zero in models A, B and C. In models Ao.oos and Aq 25 it 
was found to be around 40 per cent, while for model A0.25 
almost all the mass input is accreted. The area of the disc 
over which accretion took place in models Aq.oos and Aq 25 
was chosen so that the accretion rate and wind flux were 
nearly in balance. In general, /iacc tends to be close either 
to zero or to unity: we leave the precise elaboration of this 
dependence as a function of physical parameters for a future 
paper. 

We show the mass fluxes for models Aq.oos, A0.25 and 
A0.25 in Fig. |2^. Fig. pij a shows how the accretion region 
in model Aq.oos, which is relatively small, diverts a fraction 
of the mass which would otherwise have been driven out 
in the wind, perhaps slightly extending the length of the 
variable period compared to model A. In contrast j_if accre- 
tion occurs over a large area (Model A0.25, Fig. Bab), the 
outgoing wind is suppressed. The solution soon relaxes to a 
near-equilibrium between mass input and accretion, modu- 
lated at a 10 per cent level by the instability of the termi- 
nation shock and accretion. The evolution of model Aq 25 
(Fig. p3|c), in which accretion takes place in the outermost 
regions of the disc, is strikingly similar to that of Model 
A0.QQ5. This can be understood, as already discussed, on the 
basis of morphological arguments. 

To supply the luminosity of a QSO, L ~ 1.3 x 
10'"'/EddMh,8 erg s~^ requires accretion at a rate 



M ~ 2.3 



/e 



^acc , — 1 



Mh.gMgyr \ 



(28) 



where 77acc = 0.1r;acc,-i is the accretion efficiency of the 
black hole. The accretion rate is given, over the long term, 
as a fraction of the mass-loading rate from the stellar cluster 
by /iacc, which is determined by the global hydrodynamics 
of the nucleus. The fraction of input mass which must be 
accreted is given by 



/^acc 



0.23 



/Edd 



Jyacc-lQ- 



lOMh 



(29) 



Perry (1993a) argues that /^acc — 0.5. The accurate cal- 
culation of this fraction is one eventual aim of this work: 
some values of ^acc are given in Table ^ To the accuracy 
of the present treatment, models Aq.qqs and Aq 25 are self- 
consistent, and closest to the case considered by PD and 
Perry (1993a). 

In a Seyfert nucleus residing in a rejuvenated QSO, the 
ratio of Afh to Md may be substantially larger than 1:10. 
From equation (^^, a self-consistent model (with r/acc < 1) 
must be significantly below the Eddington limit unless gas 
is provided to the black hole from non-stellar sources or the 
stellar lifetime is very short. 

From Table ^ we see that rates for these models 
range between 3.2 and 8.6MQyr~^. Comparing these with 
the accretion rate required to fuel the central luminosity, 
4.6r;acc,-i Mq yr~^, we see that these models are reasonably 
self-consistent, i.e. the accretion rate is sufficient to power 
the central luminosity for a reasonable assumed radiative 
efficiency. If we impose the condition that the flow is self- 
consistent over the short-term, then the flow may to relax to 
one of these solutions, however sensitive the dependence of 



© 0000 RAS, MNRAS 000, 



Hydrodynamics of the ISM in Symbiotic AGN 41 



a) 



Model A 



b) 



Model B 



- Loading 
Outward advection 
Trapped outflow 



- Loading 
■ Wind 



1,5x10^ 
Time/ yr 



3x10^ 4^10^ 
Time/yr 



Figure 22. Mass processing rates of models a) A, b) B. Model A soon reaches a steady equilibrium between mass input and mass 
outflow. In Model B, mass is lost episodically from the nucleus. While as much as 60 per cent of this loss does not have escape velocity, 
we argue in the text that the effect of assuming this mass is lost to the grid will be small. Model C has no appreciable mass loss in 10^ yr. 
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Figure 23. Mass processing rates of models a) Ao.ooSi b) Ao.25, c) Aq 25- Models Aq.oos and AJ, 25 reach a steady equilibrium solution, 
very similar to Model A but with a fraction of the mass input diverted into accretion. Model Ao.25, on the other hand, has a global 
structure rather similar to Model C. Only a small fraction of the gas escapes accretion, and the central wind is consequently too weak 
to drive out through the cluster core. The mass input and accretion rates soon reach equilibrium, modulated by near-periodic instability 
of the accretion/ wind interface. 
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Figure 24. Mass within the grid for Model B, cf. Fig. ^. 
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Figure 25. Wind power in Model A. 



Atacc on flow parameters. Note, however, that delays in the 
conversion between mass accretion onto the disc and the 
output of radiation are hkely to imply that the equilibrium 
is, in fact, unstable. 



If i.11 relevant dimensionless quantities, incl uding the 



was converted into outflow power in the models presented 
here was in Model A, where the efficiency was I/wind/isoi — 
8x 10"*. Mass loss rates from the nucleus to the surrounding 
galaxy vary from effectively nil to 8.6 Mq yr~^ (or, rather, 
that the entire of the mass input to the cluster is lost). 
Rather lower values are typical in most flows. By scaling 
the models, as discussed in Section 3.2, we see that the 



wind efffciency varies as <j>\ /t : it would seem to be able 
to be increased to an arbitrary degree. However, the total 
efficiency obtainable is limited in reality by the Thomson 
optical depth of the nISM (which limits the column density, 
oc 4>\/t) and by the onset of relativistic effects as dynamical 
velocities (cx \/t) approach c. The values of the parameters 
for Model A were also chosen to reflect those likely in real 
nuclei, except that the driving radius of the flow (i.e. the 
smoothing region) is larger than is likely to obtain in real 
AGN, so values rather greater than -Lwind/ifioi ~ 10~^ may 
be possible. 

The effect on the galaxy of the winds from a symbiotic 
nucleus will be difficult to distinguish from those of winds 
driven by larger-scale galactic sta rbursts or of relativistic 
jets driven frorn smaller scales (e.g. [Veilleux et al. 1994 



Col- 



bert et al. 1996). Heckman et al. (1990) observed starburst 
superwinds with wind powers up to several lO'*'^ erg s^^ , with 
overall iwind/^Boi ~ 10~^ and wind velocities from a few 
100 km s"'^ to above 1000 km s"'^ (although again the winds 
in some of their sample may be driven, in part, by AGN). 
Whether the source is a starburst or an AGN may be distin- 
guished by whether the starburst is strong enough to drive 
the wind, and to produce the observed gas velocities, or 
whether the morphology of the wind bubble suggests an un- 
resolved, dynamically distinct source (e.g., if it is misaligned 
with galactic structure). lU-coUimated winds may be dis- 
tinguished from jets by the length-to-breadth ratios of the 
structures which result, and by evidence for relativistic mo- 
tions (e.g. synchrotron emission or superluminal motion). 

A wind with the strength we find can have a major effect 
on the ISM of the galaxy and on its surroundings. As the 
wind-blown bubble moves out, it will generate dense shells 
of gas which may be important for narrow li ne emission; 
eventually it may evacuate the ISM comple t ely (jPyson, Falle 
Perry 19 8 lt iFalle, Perry fc Dyson 198l|; traylor, Dyson & 



Edding 

L/M \mes'^!&'TJ^'^^!^^'7n^LcpenSsr^'or^ii^en^rT^(r 

t imescales cliosen (because L/ Al oc M^^/ <pMci\ by contrast, it 
L (X }? and A = r, then L/M oc r/0, and any solution with 
finite accretion can be made self-consistent). Note, however, 
that delays and instabilities in the processing of gas through 
the accretion disc may mean that, on the short timescales, 
not all AGN are strictly self-consistent, in the formal sense 
defined above. 



7.6 Wind powers 

The average wind powers found for 'mature' flows in the var- 
ious models presented here are included in Table |. As can 
be seen by comparing Fig. |25|to Fig. p2|a, the time variation 
in the power of the wind out from the nucleus was very sim- 
ilar to the variation in mass loss rate. By contrast, in model 
B, only a very little of the mass lost to the grid has escape 
velocity, so the value quoted in Table |^ is very uncertain. 
The highest efficiency with which the central luminosity 



1997 



Axon 1992; Smith, Kennel & Coroniti 1993 



iitcffon et al 



As it blows outwards, it will be an important source of 
energy and processed materia l to the intergalactic medium 
(|Voit 1994 ^ilk fc Perry 1999[). 



8 CONCLUSIONS 

The central thesis of this paper is that stars are of fun- 
damental importance in the nuclei of active galaxies. The 
symbiosis between a compact nuclear starburst stellar clus- 
ter and a massive black hole can self-consistently explain the 
properties of active nuclei. 

These two simple ingredients - a spherically symmetric 
starburst stellar cluster and an accreting black hole - con- 
spire symbiotically to produce flows of immense variety and 
complexity. Where the gravity and radiation field of black 
hole and accretion disc dominate at the very centre of the 
nucleus, a bipolar circulation is found. Outflows along the 
axis perpendicular to the accretion disc coexist with inflow 
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in a broad flattened structure above the disc. The interac- 
tion between this central flow and that in the region where 
the gravity of the stellar cluster dominates the dynamics 
results in highly complex nuclear wide flows. These result- 
ing parsec-scaled nISM bipolar flows range from episodic, 
explosive, percolating structures to well-coUimated, narrow, 
hypersonic hydrodynamic jets. The jets arise when the flow 
through most of the core is quasi-hydrostatic, while episodic 
'bubbles' occur when the Keplerian velocity of the cluster is 
substantially supersonic. It is important to note that these 
dramatic structures arise as a result of pure hydrodynamics 
in a relatively simple system; they do not require the in- 
tervention of radio-related nuclear phenomena. Our results 
dramatically illustrate why broad emission line studies have 
consistently failed to identify any simple, global flow pat- 
terns. The real flow structures in AGN are inevitably more 
complex. 

Our simulations confirm the overall flow properties of 
derived by PD and Perry (1993a). We also flnd details that 
their treatment overlooked. The gravitational attraction of 
the stellar cluster can lead to gas being retained within the 
nucleus for many dynamical times, with the line emission 
likely to be dominated by gas in a central plume which is 
retained by the gravity of t he clu ster. 

We propose (see Section 3.4.2) a classiflcation of the flow 



structures in terms of dimensionless dynamical parameters: 
the velocity dispersion of the cluster and the outflow veloc- 
ity from close to the black hole, both expressed as fractions 
of the sound speed in the nISM. These relate to the gravita- 
tional and thermal energies in the nISM. This classiflcation 
scheme, which is analogous to the Hertzsprung-Russell dia- 
gram for stars, will be broadened by the various secondary 
effects which depend on the dominant parameters. But the 
morphology of the flows are likely to remain similar and 
analytic estimates of flow velocities and densities are still 
readily calculable. 

We have discussed other observational diagnostics of 
these flows. The mass-losing stellar cluster acts as an 
optically-thin fuel reservoir which chemically enriches the 
nISM. When we consider the properties we expect in the 
broad emission line region (BELR) we find that the nISM 
itself has a multi-component structure. This fact, combined 
with the close relationship between the accretion disc and 
the kinematics of the global wind helps explain the observed 
multi-component structure of the line spectra and the ob- 
served correlations between low- and high-ionization broad 
emission lines. The varied kinematics of the flows seems to 
explain the contradictory evidence that has been found for 
infiow, outflow, chaotic or rotating flows by observers. The 
winds driven from the nuclei will also affect the surrounding 
galaxy and the IGM. 

Our study clearly shows that active galactic nuclei will 
only be fully understood once the symbiotic relationships 
between the black hole, the nuclear starburst stellar cluster, 
and the wider galaxy are considered. Our conceptually sim- 
ple Symbiotic Model provides a self-consistent explanation 
for the observed complexity of active galaxies. 
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APPENDIX A: DETAILS OF THE PHYSICAL 
ASSUMPTIONS 

In this appendix, we discuss in more detail the assumptions 
of our model. The presentation progresses from the struc- 
tures within the nucleus - the black hole, starburst stellar 
cluster and accretion disc - to the properties of the radia- 
tion field, thermal equilibrium of the nISM and treatment 
of accretion. 



Al Black Hole 

Since the bolometric luminosities of most AGN can be de- 
termined to order-of-magnitude or better, the combined is- 
sues of black hole masses and Eddington ratios are strongly 
linked. 

Indications of the masses of dark objects in nearby 
galaxies have been obtained from the brightness distribution 
and stellar kinematics in their cores. The masses inferred to 
reside in dark, unresolved cores vary from 2 x 10^ Mp in the 
Galaxy and M32 to around 3 x 10^ Mq in M87 ( [Kormendy 
fc Richstone 1995 ) . The sample of galaxies observed is at the 
moment rather sparse and heterogeneous, but the presence 
of these 'dark objects' suggests that there may be enough 
nuclear black holes of sufficient mass to satisfy our model. 

Limits on have also been derived on black hole masses 
in active AGN. Some attempts have been made to estimate 
the mass f unction using simple estimates of the black hole 



mass (e.g. Padovani, Burg & Edelson 1990), but has only 



been very recently that radio observations of maser emission 
and detailed line profile studies have begun to give well- 
constrained dyna mical masses for black holes, in the range 



10'-2 X lO'^'MiT) (GreenhiU et al. 1996; Lasota et al. 1996 



Newman et al. 1997) 



A2 The Stellar Cluster 

Dynamics We base our assumptions about the dynamics 
of the nuclear stellar cluster in the models of MCD. They 
modelled the evolution of the stellar cluster around a central 
black hole by a multi-mass energy-space Fokker-Planck code, 
including the effects of tidal disruption, physical collisions 
between stars and of stellar evolution. 

As our aim is to model the hydrodynamics of the nuclei, 
we only include their results in a schematic fashion. MCD 
find that the dynamical evolution of the stellar clusters oc- 
curs, for the most part, in the cluster cores. For initial stellar 
densities below 10^ M© pc~^, the dynamical evolution of the 
cluster is dominated by tidal disruption, leading to a core 
density law s = 7/4 (as predicted by early models of this 
process) , while for higher stellar densities collisions dominate 
the core evolution and mass loss, giving a core power-law of 
s = 1/2. 

Beyond the core, the density law in the models pre- 
sented by MCD remains roughly constant, except for the 
effects of stellar evolution. In their initial Plummer model, 
MCD take the halo power-law index to he h = 5. Beyond 
10 pc, the stellar density will be dominated by the normal 
stellar distribution of the central regions of the galaxy. The 
precise form of the halo power law and the surrounding 
galactic bulge are not well determined. However, they will 
have only a small effect on the hydrodynamic models of the 
flows inside the BELR that we present here. 

The relationship between the enclosed mass and ra- 
dius in the Galactic cen tre between 1 and 100 pc is roughly 
Af ~ 3 X Kfrl^MQ ( [Genzel fc Townes 199"^ ), implying 
stellar densities n* ~ 10''(M*/1 MQy'r'^-" pc"^. Within 
the Galactic centre, the cluster of 10* or 10^ Mq which we 
discuss here would be a dynamically independent system. 
Indeed, any such cluster would dominate the dynamics of 
at least the central 100 pc of the Galaxy. In NGC 3115, a 
dense stellar cluster of mass ~ 4 x 10^ Mq , distinct from 
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the core of the galaxy, appears to surrou nd a dark object 
(perhaps black hole) of mass ~ lO'' M© (Kormendy et al 



1996a0 although such distinct stellar clusters are not found 



around all candidate massive dark objects (Kormendy et al 
1996b|] 

HST observations have provided e vidence for stella r 
densities at least > lO'^ M© pc"^ in M32 ( [Lauer et al. 1998| ) . 
It is more difficult to observe such dense compact clusters in 
more distant galaxies, but the properties at larger radii seem 
similar to M3 2, and nuclear clusters of 10^ L© are certainly 



observed ( pL/auer et al. 1995 ), although there is no strong ev- 



iden ce f or stellar densities above 10" M© pc~'^ ( Faber et al^ 
1997 |). [ Faber et al. find that the density decreases quite 
strongly as the mass of the cluster increases, while the ve- 
locity dispersion in the clusters varies only weakly; they also 
suggest that the cluster masses in these objects are between 
3 and 6 times those of the massive dark objects they har- 
bour. A cluster with mass ~ 6 x 10^ M© contained within 
6 pc of an 4 X 10^ M© black hole is seen in the Circinus 



gala xy, the most nearby known Seyfert 2 (Maiolino et al 



199S 



The clusters which we discuss here would only be re- 
solvable in nearby galactic nuclei (although their gravita- 
tional effects would extend over a larger radius). While the 
nucleus was active, their presence would be masked by the 
active QSO and at later epochs much of their mass may 
have been lost to the nucleus in a wind or accreted by a 
black hole as a result of stellar evolution or stellar collisions. 
Clusters with these high densities are a requirement of the 
PD model if stellar mass loss envelopes are to generate the 
required covering fractions of cool gas within the BELR. 
However, clusters at least as dense are required by many 
mod els of the formation of the central massive black hole 
(e.g. Begelman & Rees 1978; Quinlan & Shapiro 1990), so 



it is not unreasonable to consider the properties of a black 
hole surrounded by the flotsam of its formation. 

The stellar cluster models of MCD matter of 

computational necessity, spherically symmetric. In the ini- 
tial hydrodynamic models presented here, we follow this as- 
sumption. In general, if the young star cluster were to form 
as a thin disc, non-axisymmetric instabilities will thicken its 
central regions substantially within a dynamical timescale, 
to generate a bulge-disc structure reminiscent of a spiral 
galaxy. If the clusters remain aspherical enough to influence 
our models, they would either consist of stars on plunging 
box- type orbits or have significant net angular momentum. 
In the first case, the aligned orbits would rapidly be ran- 



domise d by the influence of the central black hole (Gerhard 



& Binn3y 1985). In the second case, for self consistency, we 
should then also include the angular momentum of the flow. 



a topic we leave for a subsequent paper. 



Even in the core of a nuclear starbur st, the pressure, 
typically lO'^Kcm"^ ( |Suchkov et al. 1996| ), is significantly 
less than the pressure in the BELR of AGN. We have chosen 
to neglect this extranuclear gas in the simulations presented 
here, as the rate of mass input into the central parsec-scale 
region will depend strongly on the transport of angular mo- 
mentum. Detailed simulations of the gas dynamics of galax- 
ies with active nuclei will allow us to correct this in future 
papers. 



Stellar mass loss The manner in which mass is lost from 
stars and added to the flow has an effect on our simulations 
on various scales. The overall normalization of the stellar 
mass loss, within our assumptions, sets the overall normal- 
ization of the density, but does not alter the kinematics of 
the flow. Variations in mass input across the core (for in- 
stance, if the stellar cluster is oblate rather than spherical) 
will alter the form of the flow, but we argue (in Appen- 
dix ^^) that this will not have a major influence on the 
structure of the flow. The graininess of the mass input, on 
the scale of the individual interactions between stars and the 
global ISM, is determined by the mode of stellar mass loss 
(through supernovae, stellar winds, collisions, etc.). This can 
have consequences for the thermal structure of the flow. Per- 
haps the most important consequence of these interactions 
will be the production of c ool gas which will emit optical 
and UV emission lines (PD). 

In all the models presented above, we take the input rate 
of gas into the cluster to be proportional to the stellar den- 
sity. The scaling chosen, q — Qpf, with Q = 3 x 10~^® s~^ ~ 
10~*yr~^, is appropriate for mass loss from a young stel- 
lar cluster with a lower mass cut off above 4M0 (ShuU 
1983; WiUiams & Perry 1994, in particular their Fig. 3b). 
In this way, our assumptions differ from those of MCD who 
used a broad IMF, extending from stellar masses as low as 
M* =0.3 Mq . For IMFs extending to such low masses, coUi- 
sional mass loss exceeds stella r mass loss at the very smallest 
radii ( [Murphy fc Perry 199^ . 

Rieke et al. (198C ) suggested that IMFs with lower mass 
cutoffs above a few solar masses were necessary to explain 
the high supernova rates and 2 nm luminosities for the de- 
rived dynamic al mass in the bes t studied starburst galaxy . 



M82 (see also Rieke et al. 1983; Doane & Mathews 199J 



Bernlohr 199; ). Studies of the variation of the mass-to- light 
ratio and metallicity with galaxy mass for ellipticals also 
suggest that 'bimodal' st a r formation occur s in some cir- 
cumstances (Larson 1986; Zepf & Silk 1996). High-biased 
IMFs may be particularly favoured in the highly sheared 
environment of a galactic nucleus, if the increase in pre- 
collapse density required for higher shears is overweighed 
by increased internal energy resulting from turbulence and 
heating of the cloud gas (cf. Silk 1977, 1995; Larson 1985). 

It should be noted that the obs ervational e vidence for 
high-biased IMFs is not clear cut (scalo 199C). Detailed 



models of the distribution of dust can explain the broad 
band s pectra at the required mass-tq - light ratio with normal 
IMFs ( JDevereux 19891 ; jcalzetti 19971 ). In M82, the distribu- 
tion of sizes of the radio supernova remnants (and the un- 
changed population over a 10 year period) s uggest that the 



supe rnova rate was initially overestimated (Muxlow et al 



1994). While so me Galactic open clusters are deflcient in the 
1-10 M0 range (Scalo 199C), in the best nearby analogue of 



a starburst system, the central region of 30 Doradus, imag- 
ing studies flnd populations of stars do wn to SMp) with no 



sign of differences in the global IMF (Hunter et al. 199E 



Brandl et al. 1996) 



Note that we assume the formation of the cluster occurs 
on a timescale of a few million years or less. Many of the 
constraints described on the number of quiescent low-mass 
stars in the starburst IMF become yet more stringent if stars 
have formed over a period significantly longer than the main 
sequence lifetime of the massive stars, since they depend 
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not on the rate of formation but on the total mass of low- 
mass stars. There is good evidence that the formation of 
massive stars in OB associations is coeval to within 20 Myr 
( 3hull & Saken 1995), with the massive stars destroying their 
nascent star forming region (perhaps encouraging further 
generations of massive stars to form at a larger distance) . It 
is not unreasonable for closely synchronized star formation 
to have occurred in the compact stellar clusters we discuss, 
since the dynamical timescales are short, and even at 10'' K 
the sound-crossing time of the nucleus is only 10'' yr. 

Less high-biased stellar IMFs will decrease the mass in- 
put rate, and hence the density scale will be lower. Different 
modes of stellar mass loss, driven by stellar collisions, tidal 
disruption or the effects of the nuclear radiation, will boost 
the mass input rate, albeit with a different spatial distribu- 
tion. It does become difficult to produce the required lumi- 
nosity in such low-loading models, unless the cluster mass is 
extremely large. 

Stellar mass-loss through collisions is a threshold pro- 
cess: the collisions lead to disruptions when the impact ve- 
locity is greater than the surface binding energy of the stars 
(cf. MOD). If the velocity dispersion varies with radius, this 
implies that coUisional mass loading can occur from any in- 
dividual class of stars wherever the dispersion is greater than 
the stellar escape velocity. From the above discussion, we see 
that within the core, head-on collisions are likely either to 
be disruptive through most of the core (where «k > «*), or 
not at all. The coUisional mass-loading term may then be 
given by 



(Al) 



when > u*, and zero otherwise. For our purposes, this 
only need be applied to main sequence stars, since giants 
will soon lose their extended envelopes as a result of stellar 
winds and supernovae. CoUisional mass loss will thus only 
be important when Tcoii = I/tidk^* < t*. Since for a uni- 

in7/2 3/2 

form cluster Tcq w — s yf; r oughly independent 



of stellar mass (Perry & Williams 1993), collisions are only 



likely to be an important factor in stellar mass loss only for 
the small fraction of stars in the central d ensity cusp, or in 



Perry 



the den sest of nuclear clusters (cf. MCD, Murphy 
199S|)~] 

Close to the central black hole there is an exceptional 
region of increased velocity dispersion. An 'extended loss 
cone' may form, for which the vigour of the stellar environ- 
ment in the velocity cusp performs a similar destructive role 
to the high tidal field of the black hole at yet smaller radii. 
This region, together with that in which mass is input by 
tidal disruption, is not resolved by our simulations. Central 
collisions could most easily be catered for by an unresolved 
central mass addition term, to be traded off against the con- 
densation of mass onto an accretion disc. 

Effects of Loading Distribution We assume spatially 
smooth mass loading. The validity of this assumption de- 
pends on the efficiency of mixing between mass input and 
global How, and on the manner in which stars lose mass. Su- 
pernovae, tidal disruption and head-on collisions will add 
mass in very localized regions, but initially with a wide 
range of velocities. Smooth mass loading may be a rather 
rough approximation where these effects dominate the mass 
input: indeed in Seyfert nuclei, individual supernovae will 



be rare, but may be able to blow away much of the nISM 
(PWD). However, sufficiently many supernovae should occur 
in a QSO that smooth loading is a reasonable first approxi- 
mation (PD). 

Mass loss by stellar winds will result in a smoother spa- 
tial distribution of mass-loading, as many sources are ac- 
tive at any time, an d each one is active for many dynamical 
timescales. Maeder (1992) suggests that for solar-abundance 
stars over 25 M©, most of the stellar mass is ejected in be- 
fore the supernova explosion. For stars of 40 M© or over, the 
mass of the pre-supernova star is less than 8 M© . The frac- 
tion of mass lost pre-supernova may be further increased by 
the higher metallicities expected in nuclear stellar clusters. 
Glancing collisions, particularly with the envelopes of giant 
stars, will also load the region in a gradual fashion. 

The large-scale distribution of loading is important 
when rq <: pv [from comparison of terms in the equations of 
mass conservation and motion, eqs ( |l^ and ([l4[)]. Hence, the 
details of the mass input distribution are important in the 
outer accretion fiow and close to stagnation points (where 
II ~ 0). If, in the absence of accretion, the central convective 
region has size rcon at time t (where rcon saturates at rc 
at time tfiii when the wind mass loss becomes equal to the 
mass input, after which time we should substitute t = tfai in 
the following expression), we find that within the convective 
structure 



pv 



(A2) 



where fdyn = Tcon/v- Unless the onset of outfiow is prompt 
(i.e. within a dynamical timescale), the precise form of mass 
loading distribution through the cluster core will be unim- 
portant. Even if the outflow is prompt, its structure will be 
simple: the central outflow will still be independent of the 
form of the mass loading, and the outer structures will be 
ballistic (cf. Model A). We will investigate these assumptions 
in more detail in a future paper. 



A3 The Accretion Disc 

In the limit of low external pressure, a thin accretion disc 
will have a corona-wind structure similar to stellar wind 
models. The rate of mass loss in the wind would depend 
on the incident flux (whether direct from the nucleus or 
back-scattered from the nISM) and the structure of the disc 
(self-gravity may well lead to a rather clumpy structure at 
these radii), and will be influenced by magnetic fields and 
angular velocities. 

However, the models we have calculated show that the 
flow pressure tends to be high in the plane of the disc. The 
radiation intensity is low in the disc plane. Conflned, high 
pressure gas may cool, adding gas to the disc. The impor- 
tance of this process depends on the ratio of the cooling time 
of the shocked gas to the flow time through the high-pressure 
region. In one limit, essentially all the gas will 'collapse' onto 
the disc surface. In the other, mass addition from the wind 
would have a negligible effect. 

To treat these processes accurately requires detailed 
heating/cooling models of the flow in and around the disc. 
The length- and time-scales important for such models are 
difficult to inc lude in a global model for the nuclear flow. 
Woods et al. ( 1996 ) treat the evaporation of gas from the 
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disc surface, using an adaptive grid technique. In their 
model, the regions of cool gas are localized near the disc 
surface: in general, where the unstable cooling of the gas 
may be important, complicated turbulent structures would 
quickly overwhelm any adaptive grid model with irrelevant 
detail. To study these cases, local results on the detailed 
structure of the flow must be smoothed to predict the mean 
properties of gas on a global scale. 



A4 Gravity 



The gravitational forces are dominated by the central ac- 
creting black hole (point source), and the nuclear starburst 
stellar cluster. We discuss observational constraints on the 
masses of the.se components in Appendices M and , re- 
spectively. 



We do not include the gravitational field of the inter- 
stellar gas. At any time, the integrated mass in the ISM 
interior to 1 pc is far smaller than that of the stellar cluster, 
since gas will not be retained in the nuclear ISM for longer 
than the average stellar lifetime ^ lO" yr. If the mass of the 
ISM were gravitationally significant, the nucleus would be 
optically thick to Thomson scattering and X-ray absorption, 
and the ISM would rapidly cool by bremsstrahlung. 

The accretion disc is limited by its vulnerability to lo- 
cal self-gravitating instabilities, to have a density no greater 
than the tidal density. 



Mh + Mci(<r) 
(4/3)7rr3 ' 



(A3) 



and hence it can contribute no more than H/r of the gravi- 
tating mass. 

For the model stellar cluster given by equation (^), the 
cluster mass is given by 



h 



-47rp*,cr(, 



(3-s)(/i-3) 
where the core mass Mc is 
1 



h-3 



Mc 



Mc = 



-47rp* 



(A4) 



(A5) 



The gravitating mass of stars inside radius r is then 
p*{r /rc)Mc\, where 



/i-3 



P*{r/rc) 



A5 Radiation 



{i-s){r/rcr 
{h~s) 



r <rc 



r > Tc. 



(A6) 



The dynamical effects of the black hole are not limited to to 
its gravitational attraction - the high luminosity generated 
by accretion will generate a significant outward radiation 
force. Accretion on to a black hole at a rate Ma will produce 
a luminosity of 



LboI ~ rja^ccMnC 



(A7) 



where r^acc is an efficiency factor, of order 0.1. In this pa- 
per, we have taken the driving to result from optically thin 
Thomson opacity only, as it is the dominant opacity when 
the nISM is at the Compton temperature. The dynamical 



effect of this radiation can be parameterised in terms of the 
angle- averaged Eddington ratio 



(/Edd) = Lboi/L-^ 



(A8) 



between the central bolometric luminosity and the Edding- 
ton limit luminosity, cf. equation ^. 

The broadband continua of many AGN are dominated 
by a quas i -thermal 'bump' at U V wavelengths (ganders 



et al. 1989; Zheng fc Malkan 1995), which is believe d to be 



emitt e d by th e central regions of an accretion disc (Shields 
197^; |CDMP|). If this is so, the intensity of the UV ra- 



diation field will have an angular de pende nce similar to 
F{e) oc (cos 61)°, with a = 1. Netzer ( |l990| ) suggests that 
using F{9) oc cosS(l + l.ScosS) better accounts for the 
effects of limb-brightening at opt ical wavelengths (Laor & 



Netzer 198£; Pension et al. 1990) find that for NGC 4151 



intrinsic continuum anisotropy (as opposed to reddening) 
must be part of the reason for the decrement in ionizing lu- 
minosity by a factor ~ 13 between that observed and that 
inferred to excite the ENLR emission in the galaxy. While 
using a = 1.5 gives an almost identical F{9) to this, we 
can also obtain an very similar fit (at the important angles 
where the net force is outward) by rescaling the hole mass 
and Eddington ratios, keeping a = 1. For example, as a test 
case, we treated a model in which a — 1.5, with all other 
parameters as in Model B, were almost identical to those for 
Model A with a = 1 (with flow velocities about 15 per cent 
smaller, as expected). Hence we will only present results the 
simplest model here. 

We smooth the radiation field close to the centre of the 
nucleus, principally for numerical reasons. However, the ra- 
diation field may be physically smoothed in many sources, if 
the nISM is optically thick to scattering close to the nucleus. 
The anisotropic component of the radiation will mostly be 
emitted within a radius roughly determined by the black 
body emission law 



(idisc/27rasBTi 



4x1/2 



0.017LifT-24PC, 



(A9) 



where asB is the Stefan-Boltzmann constant, and rin,4 ~ 
3 for observed accre tion discs (with a broad scatter, 
Zheng & Malkan 1993). The temperature profile T{r) of a 



geometrically-thin, optically-thick disc heated either by ex- 
ternal irradiation or viscous dissipation is 



(AlO) 



It should be noted that (/Edd), /disc and Mh are not 
actually independent parameters of the hydrodynamical 
model. The net force close to the centre (taking account 
of radiative forces) can be described by an equivalent grav- 
itating mass 

Mnot{0) = Mh - [1 - /disc + 2/disc cos 6] {/Edd)Mh, (All) 

which is determined by only two parameters. However, we 
retain all three parameters, because of their physical rele- 
vance - processes outside the domain of our hydrodynamic 
models, such as the stellar dynamics close to the black hole 
or the details of the inner accretion disc do, of course, de- 
pend on the separate values of these parameters. 
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A6 The thermal state of the nISM 

The thermal equilibrium of the nISM is controlled by Comp- 
ton heating and cooling. Neglecting relativistic and occupa- 
tion number effects, the net Compton heating rate, per unit 
volume, is 



Qc 



47rncO"T 



JvhvAv — 'ik^Tc I Jv Av 



(A12) 



The ratio between the Compton cooling time and the flow 
time (~ rpc/XT7''^10" s) is 



17 



1/2 



8 X 10" 



/EddMh,8 



(A13) 



where M is the flow Mach number, and T7 its temperature 
in units of 10^ K. The Compton heating time is T/Tc times 
tc,v Only at velocities comparable to or smaller than the 
sound speed can the temperature structure of the ISM gas 
have a significant dynamical effect on the gas flows. Thus, 
adiabatic expansion can begin to cool the global gas flow 
substantially below the Compton temperature only at large 
radii 10 pc). The wind velocity will already be supersonic 
by the time gas reaches radii large enough that it can cool 
- it will coast outwards on its residual momentum until it 
hits the galaxy ISM. 

We have taken the Compton temperature as uniform 
throughout the flow. This approximation may break down 
for a variety of reasons. The variation in the contribution of 
the accretion disc compared to the direct radiation from the 
central source will almost certainly change the spectrum; 
indeed, the spectrum of the accreti on disc itself depend s 
on azimuthal angl e in some models (yun fc Malkan 1989 ), 
although not all (Czerny & Elvis 1987). It may also vary 
with distance from the nucleus as a result of, for instance, 
photoelectric absorption or oblique scattering by moving gas 
(although these processes will only be important once the 
Thomson optical depth though the ISM is close to unity). 
The temperature will also increase, at smaller distances than 
our simulation s resolve, as a result of stimu lated Compton 
scattering (cf. Rees, Netzer & Ferland 198£). 

The likely variation in Compton temperature of the gas 
with the strength of the accretion disc component may be 
investigated by varying the stren gth o f the UV bump in the 
spectrum, and solving equation (4.12) (or, in these results, 
more accurately by including relativistic correct ions to the 



cross section, as in 



Rees, Netzer & Ferland 198£). We com- 



pared the resu lts for the default spectrum given in Cloudy 
(Ferland 1996) with the same spectrum with a power- law 



interpolation under the UV bump component. Without the 
bump, the Compton temperature increases from 1.3 x 10^ K 
to 2.6 X 10^ K, and the luminosity is cut to 0.52 of its previ- 
ous value: the UV bump has a negligible effect on Compton 
heating, so Tc oc 1/(1 ~f /disc)- Thus, the variation in the 
Compton temperature with azimuthal angle is likely to be 
small enough to be negligible (if the broad-band continuum 
is isotropic). This also suggests that, taking this mean spec- 
trum to correspond to a population average, the luminosity 
of the reflected component is characteristically equal to the 
broad-band component (detailed studies of form of the X- 
ray spectrum, as det ermined bv Compton re flection, suggest 
a similar result, e.g. George fc Fabian 1991). 



Whether Comptonization is sufficient to maintain strict 
isothermality or not, it is likely that the fiow will remain hot, 
at least until bremsstrahlung cooling becomes comparable 
to Comptonization, which occurs when the fiow density in- 
creases above 



riBr 



Q 1 n5 fEddMh.a 

3 X 10 -7— cm 



I pc-1 7 



1/2 



(AM) 



i.e. the hot phase ionization parameter decreases below 
H ~ 50T7"^/^ In some of the simulations presented here, 
particularly where gas is retained for many dynamical times, 
the densities violate this limit in some places. The extension 
of the models presented here to treat the effects of gas with 
cooling no t do minated by Compton upscattering is discussed 
in Section 



7.4. 



We can, however, maintain the physical relevance of our 
models either by rescaling dimensional parameters - loading 
rates, sizes, sound speeds - or by accepting that the gas will 
cool (thereby increasing the mean opacity and the effective 
radiative driving force). We discuss these implications more 
fully in Section 7.4 - for the present, it sufficient to warn 



that scale of the densities in the models we present must be 
treated with caution: relative values and fiow topologies are 
likely to accurate nevertheless. 



APPENDIX B: EVOLUTION OF ISOTHERMAL 
FLOWS IN STELLAR CLUSTERS 

In this appendix, we derive an approximate analytic model 
for the transient behaviour of an isothermal (cs = const) 
flow in a uniform, spherical stellar cluster with no black 
hole. This model serves as a test case for our numerical mod- 
els. It also illustrates that even in spherical symmetry, the 
flows within stellar clusters can take very many dynamical 
timescales to reach equilibrium. 

We consider the case where the Keplerian velocity at 
the edge of the cluster is supersonic, vk ^ Cg. When this 
is not the case, gravitational effects are negligible, and the 
flow will relax to equilibrium on a sound- crossing timescale. 
Here, we derive how this relaxation timescale increases as 
the gravitational field increases. 

We take the stellar cluster to have a uniform core and 
no halo, i.e. s = and h — 00. The mass loading rate, q, is 
uniform throughout the cluster following its abrupt turn-on 
at time t = 0. The gravitational force is then 



— J r for r < rc 

9=\ )I,Zi2 



(Bl) 



where = GMc/tc- 

If there were free accretion at the centre of the cluster, 
the fiow would relax to a doubly-sonic equilibrium solution. 
This solution has an outward-going sonic point at a radius 

— (uK/2Cs)rc (where the flow density is p^), and an 
inner sonic point within the cluster. The ratio of mass fluxes 
between outer and inner sonic radii is 



47r (r^T) 2 p-Cs 



(B2) 



While the area ratio increases between the sonic radii, the 
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Figure Bl. A section tlirougli Model 0, shown in Fig. ^ (in which mx = 7). The points plotted show a) density b) total velocity, 
along the line r = 2. In a), dashed lines show the density derived from the hydrostatic model for the core and the ballistic model for 
the accreting halo. The velocity profile in the accreting halo model is shown in b). The numerical results agree well with the core model 
inside the shock at r ~ 0.5 pc, and with the accreting halo model beyond it. Small deviations in the core probably result from vibrations 
or some supra-thermal support by turbulence. 



density decreases exponentially in the near-hydrostatic re- 
gion between the sonic points and the total mass flux ratio 
is small, so little mass is lost to the cluster at this stage. 
The stagnation radius, ro, (where the flow velocity is zero) 
is only just inside the edge of the cluster core. 

To determine the flow within the cluster, we use the 
mass and momentum equations, equation ([l^ (integrated, 
assuming the stagnation radius is coincident with the edge 
of the cluster) and (|l4|), to obtain 

ir'^pv = -q{rl~r^) (B3) 
Av vie q 



5-r« iv -fcJ+2 . (B4) 

p ^ ' r 

Below, we use dimensionless variables m — v/cs and x — 
r/rc (so that rriK ~ vk/cs is the Mach number at the Kep- 
lerian velocity). 

We approximate the accretion flow by free-fall collapse. 
This will be strictly applicable only when there is an inner 
sonic point, which we flnd is at Xs given by 



(m|xs - 2)(1 - x^) 



6x^ 



(B5) 



This equation has no solutions for Xs between and 1 if 
rriK < 2 X 2^''^ ~ 3.17. For niK > 3.2, the relevant solu- 
tion, x~ , may be approximated by expanding equation (B5) 
about x~ = 1: 



(B6) 



Assuming x < x^ 
to obtain 

dm 2 Sx'^m^ 

m — = -m-i^x + 

ax 1 — x'^ 

Taking m 



1, we take Cs — in equation (B4) 



(B7) 



at s = 1 we flnd the inward velocity is 



^ l^ 3^-l / 9- 20x2 + 16x5 -5x8 
v = -v^x{l-x) d — . (B8) 

When the accretion rate at the centre is less than that 



at which mass is added throughout the flow, gas will pile 
up in a hydrostatic core, surrounded by an accretion shock. 
The density distribution in the core (assuming u ~ here) 
is 



p = ps exp 



2?j 



(B9) 



where ps, Xs are the values just inside at the accretion shock. 
These may be found, at any time, by equating psC^ with the 
ram pressure of the accretion flow, pv^, and by setting the 
mass accumulated in the core equal to the total mass input 
prior to t: 



An 3 



p{x)4tvx dx. 



(BIO) 



The fit between equations (B8) and 



|) and the time- 
dependent numerical model. Model above, is illustrated 
in Fig. ^ 

Together, these equations (Bt and Bl(j) give the time 
at which the accretion shock reaches Xs, t( a function 

of rriK- The function t{xs) is zero at Xs = and Xs — 1, with 
a maximum at t = tmcrgo at a value of Xs close to the edge of 
the stellar cluster for large rriK. At times later than tmorgc, 
the inner sonic point merges with the accretion shock, and 
supersonic inflow ceases. After this time, the stagnation ra- 
dius begins to move inwards within the cluster (when equi- 
librium is flnally reached, the stagnation radius is at the 
centre of the cluster). 

For large ttik, the timescale over which the hydrostatic 
core flUs the central stellar cluster can be approximated by 



(Bll) 



where Ik = Tc/vk is the free-fall timescale within the stel- 
lar cluster. When rriK = 3, the filling timescale is roughly 
25 times the free-fall timescale, rising to roughly 350 times 
for rriK = 4 and very rapidly thereafter. A yet longer time, 
roughly AtK/'m'^e'"^'^ , is then taken to fill the region be- 
yond the cluster core out to the outer sonic radius, rcrn^/2. 
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Only after this time, far longer than any of the obvious dy- 
namical timescales, is equilibrium finally achieved through- 
out the subsonic region. The wind through the sonic ra- 
dius, which has gradually been strengthening, then removes 
a mass flux equal to the mass input in the core. 

This paper has been produced using the Royal Astronomical 
Society/Blackwell Science style file. 
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